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ABSTRACT 

In  recent  years  the  Box-Jenkins  method  has  become  a popular 
technique  for  forecasting  future  behavior  of  a time  series.  Cnee 
the  forecast  model  is  known  the  method  is  very  easy  to  employ  and 
adequate  computer  packages  are  available  for  most  purposes.  Un- 
fortunately the  problem  of  determining  the  appropriate  forecast 
model  has,  for  models  of  any  complexity,  been  one  of  the  major 
stumbling  blocks  to  the  user  of  this  method. 

In  this  paper  a satisfactory  solution  to  that  problem  is  ob- 
tained and  it  is  demonstrated  by  numerous  examples  how  this  greatly 
enlarges  the  class  of  data  sets  which  can  be  adequately  modeled  by 
autoregressive-moving  average  models.  This  new  approach  is  suffi- 
ciently unequivocal  that  most  users  will  find  it  easy  to  implement. 


^This  research  was  sponsored  by  ONR  Contract  #N0C014-75-C-0439. 
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I . INTRODUCTION 

In  [7]  H.L.  Gray,  A.G.  Houston  and  F.W.  Morgan  have  studiad  in 
some  detail  a new  spectral  estimator  referred  to  as  the  G-Spectral 
Estimator.  In  that  paper  they  pointed  out  that  the  prohlem  of  se- 
lecting the  appropriate  G^  - transform  for  spectral  estimation  is 
almost  equivalent  to  the  problem  of  selecting  p and  q in  modeling 
an  autoregressive-moving  average  process  of  order  (p,  q) , i.e.  an 
ARMA  (p,  q)  process.  Moreover,  it  was  also  suggested  and  implicitly 
demonstrated  that  the  "R  and  s arrays"  (defined  later)  which  are 
used  to  calculate  the  G-Spectral  estimator  could  be  utilized  to  de- 
termine p,  d and  q in  an  ARIHA  (p,d,q)  process.  In  this  paper  that 
suggestion  is  explored  in  depth  and  a new  set  of  criteria  which 
uniquely  determines  p,  d and  q and  which  is  readily  observable  is 
established.  Many  examples  are  given  and  it  is  demonstrated  that 
this  new  approach  not  only  simplifies  the  problem  of  ARMA  model 
fitting,  but  greatly  extends  the  range  of  models  that  can  reasonably 
be  identified.  This  is  particularly  true  in  the  case  of  the  non- 
stationary model  where  the  method  is  utilized  to  identify  and  re- 
move nonstationarities  due  to  real  or  complex  roots  on  the  unit 
circle. 

In  the  following  section  we  set  down  the  definitions  and  theorems 
necessary  to  establish  the  results  referred  to  above. 

II.  DEFINITIONS  AND  THEOREMS 

Definition  1. 

Let  m be  an  integer,  h > 0,  and  f be  a real  valued  function. 

Also  let 
f - f(mh), 
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It  has  been  shown  by  Pye  and  Atchison  (10)  that  R (f  ) and 

n in 

S^(f^)  can  be  calculated  recursively  by  the  following  relations. 
Define 


SQ(fjj|)  = 1,  m = 0,  +1,  +2,  ... 

= V ” 0'  ^2 


Then 


n+1  m 


n \ r ^ m+l^  , , 

'^n^^m+l^^  S (f  ) " 

n m 


(4) 


and 


n m 


(5) 


for  n » 1,2,  ...  and  in=0,  +1,  +2,  ...  . In  the  future,  when 

they  are  defined  we  use  (4)  and  (5)  to  define  R ,,(f  ) and  S (f  ) 

n+l  m n m 

even  though  the  ratio's  in  (2)  and  (3)  may  be  undefined. 

In  (4)  and  (5)  we  have  tacitly  assumed  that  S (f  ) / 0 and 

n m 

R (f  ) ;<  0.  If  R (f  , , ) = R (f  ) = 0 we  leave  S (f  ) undefined, 
n m n m+ 1 n m n m 

However,  if  R (f  ) = 0 and  R (f  ^ 0 we  define  S (f  ) = + ®.  A 
n m n m+ 1 n m 

similar  definition  will  be  employed  for  R , (f  ) . 

n+i  m 

Definition  2. 

A function  f will  be  said  to  be  an  element  of  L(n,A)  over  a 
set  of  integers  I = {mQ,mQ+l, . . . } if  there  exists  a smallest  inte- 
ger n > 0 and  a set  of  such  that  f is  a solution  of  the  dif- 

ference equation 
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- 4 y 
n''m-n 


0 for  mcl. 


For  future  reference  we  now  Include  the  definition  of  an  AKMA  (p,  q) 
process. 

Definition  3. 

A stochastic  process  {X^},  t » 0,  + 1.  + 2,  . . . is  said  to  be 
autoregressive  of  order  p and  moving  average  of  order  q if 

p q 

X - E (>..x^  . + - E e._Z^  (6) 


^ W-k  * \ 

k»l 


^ ®Jc^t-k' 
k-1 


where  the  and  0j^  are  constants  and  is  white  noise.  By  uti- 
lizing the  operator  B defined  by  BX^  ■ ^t-1  '^111  usually  write 
equation  (6)  in  the  form 

0 (B)Z^ 


where 


i(i(B)X 
«(B)  » 1 - 


t 
B - 


(7) 


0(B) 


1 - 0^B 


e^B 


0 B*’. 


The  following  theorem  forms  a basis  for  the  development  of 
several  of  the  results  to  follow.  Its  proof  can  be  found  in  [7] . 
Theorem  1. 

Let  n > 0 and  suppose  S (f  ) and  R (f  ) in  (2)  and  (3)  are 

n m n m 

defined,  i.e.  the  denominators  in  (2)  and  (3)  are  not  zero,  and 


S (f  ) ^ 0.  Then  f G L(n,A)  for  m > m + n if  and  only  if  S (f  ) 
n rn  m — u n m 

is  constant,  as  a function  of  m,  for  n>  > 

III.  ARMA  (p,  q)  Processes  - The  Stationary  Case 


I The  current  most  popular  approach  to  modeling  an  ARMA  (p,  q) 

I process  is  a method  referred  to  as  the  Box  - Jenkins  method. 

I It  consists  of  initially  examining  graphs  of  the  estimated 

[ autocorrelation  and  partial  autocorrelation  to  obtain  possible 

( values  of  p and  q and  then  investigating  these  contending  models 

more  closely  to  determine  an  appropriate  model,  see  [5).  There  are 
a number  of  difficulties  with  the  method,  however  most  of  them  stem 
from  the  fact  that  even  if  the  graphs  of  the  autocorrelation,  p(t), 
and  partial  autocorrelation,  were  error  free  the  simple  in- 

spection of  their  graphs  would  not  in  general  yield  unique  values 

! 

1 


i 

I 
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of  p and  q when  p>0,q>0>  i.e.  in  the  mixed  model.  When  this 
fact  is  coupled  with  the  ambiquity  induced  by  the  necessity  of  uti- 

A A 

lizing  estimators  p(t)  and  in  place  of  the  true  values  the 
method  becomes  richly  infused  with  art.  In  this  section  we  will 
show  how  Theorem  1 can  be  expanded  in  such  a way  as  to  establish 
new  criteria  which  uniquely  determine  p and  q and  which  are  readily 
observable.  The  following  theorems  are  the  basis  for  these  remarks. 
Theorem  2. 


Let  {X^}  be  a stationary  AIUIA  (p,  q)  process  with  autocorre- 
lation p (t)  . 

Then  for  f 


* 0 

n m 
^ 0 


Suppose  that  S (f  ) and  R (f  ) are  defined  p > 0 and 
n m n m 

p , some  integer  m.,  and  some  constant 
In  0 


s„(f  ) - c, 

n m 1 

s„  (f  ,)  C 
n m^- 1 1 


m > Oq 


if  and  only  if  n 
Moreover 


p emd  m 


q - p + 1. 


C,  - (-1)*^[1  - Z ♦.], 
^ k-1 


(8) 


(9) 


where  the  are  defined  in  (7). 

Proof.  Since  {x^}  is  ARMA  (p,  q) , p(t)  e L(p,A)  for  T>q+1- 

q - p + 1 + p but  not  for  t > q and  hence  the  first  part  of  the 

theorem  follows  from  Theorem  1.  From  the  first  part  of  this  theorem 

and  the  definition  of  S (p  ) , equation  (9)  then  follows  by  simple 

n m 

column  additions. 

Theorem  3. 


Under  the  conditions  of  Theorem  2 

C, 


n m 


S (p 
n m^+l 


2 

) / C, 


m < m^^ 


(10) 


if  and  only  if  n = p and  m^^  » 
Moreover 


"J. 


(11) 


Proof.  See  Appendix  1. 
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Although  Theorems  2 and  3 are  extremely  useful,  the  following 
theorems  will  probably  be  even  more  helpful  in  the  data  that  will  be 
considered  later. 

Theorem  4. 

Let  the  conditions  of  Theorem  2 be  satisfied.  Then  for  f ■ 

m 

(-1)  p , some  integer  m.,  and  some  constant  C,  f'  0 
m 0 3 


n n 


m > 


n 3 

if  and  only  if  n ■ p and  m^  ■ q - p ♦ 1. 
Moreover 

p ^ k 

C-  - (-D^^d  - Z (-1)  ♦.  ) 

k-l 


(12) 


(13) 


Proof.  The  proof  follows  in  the  same  manner  as  Theorem  2 by  noting 

that  f » (-D^p  satisfies  the  difference  equation 
m m 

P )t  k 

(1  - I (-1)  <p.B  }f  * 0,  when  m > q 1. 

. . iC  in  ** 

k*l 

Theorem  5. 

Under  the  conditions  of  Theorem  4 


S„((-l)  P„)  - C. 
n m 4 

m.+l 

if  and  only  if  n =•  p and  m^^  * - q - p. 
Moreover 

(-l)P^^C, 


Proof.  The  proof  follows  in  the  s2une  manner  as  Theorem  3. 
Corollary  1.  Under  the  assumptions  of  Theorem  2 


. . - . ni 
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- - q - p. 

Proof . The  result  follows  at  once  from  the  above  theorems  and  the 
recursion  rule.  Although  Corollary  1 characterizes  p and  q theo- 
retically as  well  as  our  previous  results,  we  will  seldom  use  it  to 
determine  p since  its  column  behavior  is  not  as  distinct  when 
errors  are  present  as  that  described  by  our  previous  results.  ‘Ihus, 
p will  normally  be  identified  when  errors  are  present  by  talcing  p 
as  the  first  column  having  the  behavior  described  by  Theorems  2 and 
3.  In  contrast  we  shall  see  that  the  behavior  of  the  R array  is 
usually  adequate  for  determining  q.  These  remarks  will  all  be 
clear  shortly.  In  order  to  demonstrate  the  manner  in  which  the 
above  results  can  be  utilized  to  identify  an  ARMA  (p,  q)  process 
we  will  now  consider  several  examples  making  use  of  the  true  auto- 
correlation. This  is  not  to  imply  that  the  true  autocorrelation 
is  ever  known  or  used  in  an  actual  problem  but  simply  to  demonstrate 
the  true  state  of  nature  as  pertains  to  the  above  results. 

In  each  of  the  examples  which  are  to  follow  we  will  display 
arrays  of  numbers  referred  to  as  R and  S arrays.  They  are  computed 

via  (4)  and  (5)  and  arranged  as  in  Tables  I and  IL  Since  f will 

in 

be  given  in  each  case  we  shorten  the  notation  within  the  tabular 

values  to  R (f  ) = R (m)  and  S (f  ) = S (m) . 

n m n n m n 


Sj^C-l+l)  ...  Sj^(-1+1) 

S^(-l)  ... 

Sj^(O)  ...  S|^(0) 

S^(l)  ...  Sj^(l) 


s^(j)  ...  s^(j) 


TABLE  II 
1 2 


Rj^(-4) 

RjC-i) 

R^(-t+l) 

R2(-i+l)  ... 

R^(-l) 

Rj^(-l) 

Rj(0) 

Rj^(O) 

(1) 

Rj^(l) 

Rl(j) 

R^Cj) 
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TABLE  ZV 

Cl 


m 

n 

1 

2 

.. . p+l 

-t 

p(-t) 

RjC-l) 

• • * 0 

-i+i 

p(-l+l) 

R2(-1'*-1) 

■ • • 0 

• 

-q-p-l 

0 (-q-p-l) 

R2 (~q”P-l) 

• 

• 

• * . 0 

-q-p 

P (-q-p) 

• • • Nonzero 

q-p 

P (q-p) 

Rj (q-p) 

. • , Nonzero 

q-p+1 

0 (q-p+1) 

Rj (-q-p+1) 

...  0 

j 

P (j) 

Rjtj) 

...  0 

In  Tables  III  and  IV  we  have  let  f * o • The  arrays  do  not 

mm 

differ  in  their  pattern  if  f > (-1)  o • Of  course,  in  that  case 

mm  _ 

m 

Column  1 in  the  R-array  is  made  up  of  values  of  (-1)  p . In  Table 

m 

III  we  have  included  Column  p+1  even  though  it  contains  several  un- 
defined elements.  Such  elements  are  indicated  by  a factor  of  the 
form  (5“11 • That  is,  in  the  presence  of  noise  we  would  expect  the 
column  having  the  characteristics  of  Column  p to  be  followed  by  a 
highly  variable  column  since  in  practice  the  g quantity  will  be  re- 
placed by  the  quotient  of  two  small  numbers.  A similar  consider- 
ation of  the  p+2  column  in  the  R-array  suggests  the  zero  behavior 
should  essentially  continue.  Later  ex^unples  will  clarify  this. 

Also  in  later  examples  we  will  star  the  array  elements  S (-m) , 

m 

S (-m+1)  and  R (-m). 
m m+i 


( 


J 


This  is  because  the  constancy  behavior  in  the  S-array  must 
occur  symetrically  about  one  of  these  pairs  of  points  and  the  zero 
behavior  in  the  R-array  must  occur  symetrically  about  one  of  the 
R (-m)  elements. 

Example  1. 

To  demonstrate  the  previous  theorems  consider  the  R and  S 
arrays  for  the  process 

- 1.32X^_j^  + •®®’‘t-2  “ ^t  “ *®^t-l' 


Tables  V and  VI  show  these  arrays  when  p is  known  using 

CQ  ^ 

f ■ (-1)  p and  f ■ p respectively.  It  is  well  worthwhile  to 
m n in  in 

inspect  those  tables  closely  now  since  all  of  the  previous  results 
are  well  illustrated  there  emd  the  remaining  part  of  the  paper  is 
predicated  on  those  results. 

Tables  VII  and  VIII  show  the  R and  S arrays  when  p is  replaced 

. m 

by  the  sample  autocorrelation  p , i.e.  p is  replaced  by 

in  in 

T-m 


.E,  (X.-X)(X  -X) 

p - , (15) 

m T - 

.Z,  (X.-X) 
i^l  i 

where  the  X^  are  data  simulated  from  (14)  with  the  2^  N(0,1). 

In  each  case  the  arrays  are  shown  to  at  least  one  column  more 
than  necessary  and  several  rows  more  than  necessary.  In  practice 
some  preassigned  number  of  rows  and  columns  must  be  selected.  Se- 
lecting this  much  larger  than  one  expects  is  necessary  induces  no 
problem.  For  Tables  VII  and  VIII  a sample  of  size  200  was  used  to 

obtain  p . One  should  realize  that  in  calculating  the  elements  of 
m 

the  R and  S arrays  the  number  of  lags  utilized  in  these  calculations 
increases  as  you  move  to  the  right  of  the  table  and  as  you  move 
vertically  away  from  the  S^(-i+l),  S^(-i)  and  R^^j^(-i)  elements. 

For  this  reason  the  values  more  removed  from  these  points  will  be 
more  variable.  Fortunately  the  values  in  the  proximity  of  these 
points  determine  p and  q.  In  Appendix  2 a table  is  given  showing 
the  smallest  and  largest  lag  (i.e.  value  of  m)  utilized  in  calcu- 
lating individual  elements  of  the  R and  S arrays. 
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In  observing  arrays  such  as  presented  in  Tables  VII  and  VIII 
several  approaches  can  be  taken  to  identify  the  pattern  present. 

To  these  authors  the  best  approach  is  to  first  inspect  the  S-array 

Choose  the  one  which  appears  most 


p and  f 
m m 


(-D^P 


with  f 

m n m m 

distinctive  and  identify  p as  the  first  column  having  the  correct 

constant  behavior  followed  by  a highly  variable  column.  This  is 

veiry  clear  in  Table  VII  where  obviously  p ~ 2.  From  that  same 

table  it  is  also  clear  that  q-p+l«Oso  that  q ■ 1.  Notice 

that  there  are  several  different  places  where  one  could  identify 

p and  q.  For  example  the  fact  p - 2,  q - 1 is  obvious  at  a glance 

since  the  two  starred  values  in  S2(m)  are  clearly  distinct.  Also 

note  this  choice  of  p auid  q is  confirmed  by  the  fact  that 

1s2(-4)|  = 44.4542,  the  largest  magnitude  in  the  t^d3le,  and 

S^C-l)  = - 3.5896.  This  is  of  course  significant  since  if  p » 2, 


q = 1,  the  true  value  of  |S2(-4) 


“ while  S2(-l)  “ “ 


- S^CO) 


m 


TABLE  V 
(-1)  P 


R-Array 


S -Array 


Rj^  (m) 

Rj  (m) 

Rj  (m) 

m 

Sj^(m) 

SjCm) 

Sj  (m) 

.13829 

.0378 

.000000 

-10 

-2.2806 

4.4117 

-.17710 

-.0553 

.000000 

-9 

-1.7928 

4.4117 

• 

.14042 

.1693 

.000000 

-8 

-1.0863 

4.4117 

-.01212 

.2224 

.000000 

-7 

14.0844 

4.4117 

-.18295 

-.1005 

.000000 

-6 

-3.0386 

4.4117 

ut 

.37298 

.0992 

.000000 

-5 

-2.2198 

4.4117 

u 

-.45498 

-.1530 

.000000 

-4 

-1.7356 

4.4117 

00 

.33468 

.5636 

.143284 

-3 

-.9420 

4.4117 

6.2511* 

.01940 

.4758 

.299323* 

-2 

-28.3076 

5.7403* 

-1.5639* 

-.52985 

-.4701* 

.074883 

-1 

-2.8873* 

2.0857* 

-3.0000 

1.00000* 

.1708 

.000000 

0 

-1.5298* 

3.0000 

u 

-.52985 

-.3235 

.000000 

1 

-1.0366 

3.0000 

u 

.01940 

-.3832 

.000000 

2 

16.2492 

3.0000 

u 

.33468 

.1040 

. 000000 

3 

-2.3594 

3.0000 

• 

-.45498 

-.0674 

.000000 

4 

-1.8197 

3.0000 

• 

.37298 

.0683 

.000000 

5 

-1.4905 

3.0000 

• 

-.18295 

-.1512 

.000000 

6 

-.9337 

3.0000 

-.01212 

-.1151 

.000000 

7 

-12.5774 

3.0000 

.14042 

.0376 

. 000000 

8 

-2.2612 

3.0000 

-.17710 

-.0257 

.000000 

9 

-1.7808 

3.0000 

u = underfined 
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(m) 

f 

m 

R-Ari a/ 

R2  (m)  Rj  (ra) 

” Pm 
m 

m 

S^(m) 

S-Array 

S^Cm) 

SjCm) 

.1382 

-.3078 

-16 

.2806 

r3!ra'4 

.1771 

.4789 

.0000 

-9 

-.2071 

.5294 

• 

.1404 

.2014 

.0000 

-8 

-.9136 

.5294 

• 

.0121 

.1947 

.0000 

-7 

-16.0844 

.5294 

• 

-.1829 

.2940 

.0000 

-6 

1.0386 

.5294 

u 

-.3729 

1.0021 

.0000 

-5 

.2198 

.5294 

u 

-.4549 

-1.0045 

.0000 

-4 

-.2643 

.5294 

00 

-.3346 

-.5018 

-1.194030 

-3 

-1.0579 

.5294 

-3.203096* 

.0194 

-.5119 

.973987* 

-2’ 

26.3076 

1.7640* 

-.801331* 

.5298 

-1.5298 

-.243666 

-1 

.8873* 

.6409* 

- J60000 

1.0000 

.5558 

.0000 

0 

-.4701* 

.3600 

u 

.5298 

.3481 

.0000 

1 

-.9633 

.3600 

u 

.0194 

.3412 

.0000 

2 

-18.2492 

.3600 

u 

-.3346 

.6831 

.0000 

3 

.3594 

.3600 

u 

-.4549 

-.6814 

.0000 

4 

-.1802 

.3600 

• 

-.3729 

-.1999 

.0000 

5 

-.5094 

.3600 

• 

-.1829 

-.1324 

.0000 

6 

-1.0662 

.3600 

• 

.0121 

-.1369 

.0000 

7 

10.5774 

.3600 

.1404 

-.3256 

.0000 

8 

.2612 

.3600 

.1771 

.2093 

.0000 

9 

-.2191 

.3600 
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Thus  when  one  puts  together  the  total  picture  implied  by  Theorems 

2-5,  no  other  choice  of  p and  q seems  even  remotely  possible.  The 

validity  of  this  latter  statement  will  be  even  more  clear  later. 

It  should  be  recalled  that  p appears  in  Column  1 of  the  R- 

array  when  f = p and  (-1)  p appears  in  that  column  when  f ■ 

, m m m m 

(-1)  p . Thus  in  this  example  errors  in  the  autocorrelation 
n 

estimates  can  be  seen  by  con^iaring  Tables  5 or  6 with  Tables  7 or  8 . 

In  this  example  both  p cmd  q are  easily  seen  from  the  R or  S 

array  at  f = (-l)'"p  . As  we  shall  see  this  is  not  always  the 
m m 

case  and  in  fact  the  S-array  is  usually  better  for  determining  p 

while  the  R-array  is  often  better  for  q.  Moreover,  although  ^ 

is  not  always  the  case,  f = p usually  yields  better  identifi- 
— mm 

cation  for  very  high  frequency  data  while  f = (-1)  p is  usually 

m m 

better  for  low  frequency  data.  Since  most  of  the  data  we  examine 

will  be  of  the  low  frequency  nature,  in  most  of  our  remaining 

exeimples  we  will  consider  R and  S arrays  only  at  f = (-1)  p . 

in  m 

Some  remarks  should  be  made  concerning  the  R-array  before 
leaving  this  example.  Knowing  from  the  S-array  that  p = 2 we 
direct  our  attention  to  R^ (m)  column  of  Table  VII.  Since  p = 2, 
the  number  of  moving  average  terms  is  the  number  of  nonzeros 
above  or  below  the  starred  value  in  that  column.  Since  R^C-S)  is 
5 to  10  times  larger  than  the  values  near  and  above  it  we  have  at 
once  q = 1.  Note  that  the  assumption  that  R^(-2)  / 0 implies 
R^C-l)  ¥ 0 and  vice  versa.  Finally  we  should  stress  that  the 
"zero"  and  "constant  behavior"  must  take  place  about  the  starred 
values.  Such  behavior  anywhere  else  in  the  arrays  is  not  indi- 
cative of  an  ARMA  (p,  q)  process. 

Example  2. 

Consider  the  series  referred  to  in  Box  and  Jenkins  as  Series 

D and  identified  there  as  an  ARMA  (1,  0).  The  series  consists  of 

310  points  and  its  estimated  autocorrelations  are  obvious  from 

the  first  column  of  the  R-array  in  Table  IX  since  R,  (m)  = (-D^p  . 

1 m 

Clearly  from  either  array  the  process  is  a (1,  0).  Notice 
how  closely  these  arrays  follow  the  theoretical  pattern  described 


J 
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TABLE  IX 
Series  D 


R-Array  S-Array 


(m) 

R2  (m) 

R3  (m) 

m 

S3(m) 

S2(m) 

S3(m) 

-.2703 

-.0022 

.0040 

-9 

-2.1470 

-1.7608 

-.0452 

.3101 

-.0041 

.0102 

-8 

-2.1315 

-.0295 

-1.8241 

-.3509 

-.0041 

-.4220 

-7 

-2.1563 

.0429 

-4.5032 

.4057 

-.0041 

-.0004 

-6 

-2.1343 

4.5076 

-50.2349 

-.4602 

.0044 

.0046 

-5 

-2.1532 

4.0960 

-4.1915 

.5307 

-.0039 

.0337 

-4 

-2.1712 

-.6560 

5.0495 

-.6215 

-.0051 

.0801 

-3 

-2.1849 

3.6692 

60.0155* 

.7365 

.0035 

.1355* 

-2 

-2.1697 

86.8522* 

-1.8438* 

-.8615 

-.1385* 

.0042 

-1 

-2.1608* 

) .9022* 

-.2479 

1.0000* 

.0031 

.0039 

0 

-1.8615* 

4.5121 

-3.4356 

-.8615 

-.0043 

-.0276 

1 

-1.8549 

.4281 

-3.9221 

.7365 

-.0033 

-.0004 

2 

-1.8440 

3.9733 

-39.8316 

-.6215 

.0038 

.0036 

3 

-1.8539 

3.6023 

-.0260 

.5307 

-.0035 

.0060 

4 

-1.8672 

-.0383 

-.5124 

-.4602 

-.0036 

-.1154 

5 

-1.8816 

.0254 

-.8472 

.4057 

-.0036 

.0001 

6 

-1.8648 

.8468 

-174.3048 

-.3509 

-.0020 

-.0125 

7 

-1.8838 

.8213 

-7.2421 

.3101 

-.0011 

-.0035 

8 

-1.8718 

10.1191 

-6.7227 

-.2703 

.0048 

.0053 

9 

-1.8795 

2.6900 

-5.1405 

by  our  previous  theorems.  For  example,  the  "-Cj^  and  “ behavior" 
in  $2 (m)  are  as  distinctive  as  the  constant  and  zero  behavior. 
Again  the  increased  variability  following  the  stable  column  is 
apparent.  Thus  again  when  one  considers  the  entire  picture,  a 
(1,  0)  is  the  only  possible  choice.  Of  course  this  particular 
example  is  also  clear  by  the  Box-Jenkins  method.  A final  remark 


should  be  made  here. 


and 


(16) 

(17) 


If  we  estimated  C,  by  S (p  ^, ) and  C.  by  S 
1 P q-P+1  2 i 


then  the 


resulting  ij)  is  easily  shown  to  be  the  Yule-Walker  estimate  of  iji 
P . P 

given  q.  Similarly  if  C.  is  estimated  bv  S ((-1)^  ^ "‘‘P  .) 

j p g-p+1 


and  by  S ( (-1)  ^p 


) the  result  is  again  the  conditional 


I 


(i.e.  given  q)  Yule-Walker  estimate  of  Thus  in  this  example, 
since  we  simultaneously  know  p = 1,  q = 0,  we  also  have  the  Yule- 
Walker  estimate  of  given  q,  i.e. 


* _ 1.8615 

’l  “ 2.1608 


.86. 


These  general  observations  will  be  of  value  in  Section  IV. 

Example  3. 

In  Table  X the  S-Array  and  Columns  1 and  4 of  the  R-Array  are 
displayed  for  the  process 


\ - -^Vi  - ^ 

They  were  calculated  using  which  was  obtained  from  300  simulated 
values  from  (18)  with  i.i.d.  N{0,  1). 


TABLE  X 


R & 

S ARRAYS 

at  f = 
m 

(-l)''''p  for 
m 

Data 

Generated 

from 

(18) 

Sj^(m) 

S,  (m) 

S^Cm) 

4 

m 

(m) 

• • • 

R.  (m) 

4 

-.3166 

-.0730 

.1511 

-.0843 

-9 

.7040 

• « • 

.3956 

.4645 

-.0196 

.1389 

-.1621 

-8 

.4810 

• • • 

.1556 

-.1948 

-.1069 

.1537 

-.5727 

-7 

.7045 

-.0085 

.1833 

.2332 

.1532 

-.1693 

-6 

.5672 

.0234 

.0580 

-.0120 

.1530 

90.8817 

-5 

.6712 

-.0024 

-.1301 

-.1479 

.1526 

.3836* 

-4 

.7102 

-1.4827 

.3996 

-.0554 

.2829* 

.0794* 

-3 

.6178 

-3.4936' 

-.3326 

-.1075* 

-.0658* 

.1375 

-2 

.8647 

.7232 

.7327* 

.0857* 

-.1379 

1.4015 

-1 

.5771 

.0022 

-.4228* 

.0407 

-.1377 

.1837 

0 

1.0000* 

-.0205 

.4984 

.0879 

-.1346 

2.4277 

1 

.5771 

.0075 

-.2855 

.0046 

-.1345 

-.1535 

2 

.8647 

-.1279 

.1496 

.1240 

-.1264 

-.0368 

3 

.6178 

-.2833 

-.0548 

. 3387 

-.0994 

-.0157 

4 

.7102 

-.3881 

-.1549 

.0291 

-.1115 

.0334 

5 

.6712 

-.4429 

Again  it  is  very  clear  from  the  S-ARRAY  that  p = 3 and  q =•  1.  Note 
the  behavior  of  both  Columns  3 and  4 are  exactly  as  would  be  expected 
for  a ARMA  (3,  1).  That  is,  Column  3 has  the  characteristic 
"constant  behavior".  Column  4 is  highly  variable  with  S^(-2)  = Sj(-l) 
while  ls^(-5)l  is  by  far  the  biggest  number  in  the  table.  Column  4 
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of  the  R-array  also  clearly  gives  q ■■  1.  Thus  again  it  is  clear 
that  no  other  values  of  p and  q seem  even  remotely  possible. 
Example  4.  (Sunspot  Data) 

In  this  example  we  consider  the  sunspot  data  considered  in 

Box  and  Jenkins  and  denoted  as  Series  E.  The  data  consists  of 

100  data  points  referred  to  as  sunspot  numbers  and  can  be  found 

in  [ S ] . Table  XI  gives  the  sample  acf  for  this  series  while 

Table  XII  shows  the  corresponding  S-Array  at  f - (-1)  p . The 
• m m 

R-array  is  not  given  here  since  the  S-array  is  so  distinctive. 
Both  the  behavior  of  the  second  eind  third  column  are  distinctive 
here  and  clearly  we  have  p » 2,  q = 0 which  is  the  same  as  found 
in  [ 5 ] . It  is  worth  noting  here  that  Column  2 is  surprisingly 
stable  for  such  a short  series  as  this.  Possible  explanations 
for  this  behavior  are  hidden  periodicities  and  nonstationarity , 
and  are  currently  being  considered  by  the  authors. 


TABLE  XI 

P for  Sunspot  Data 
m 


m 

1-7 

1 

.806 

.428 

.070 

-.169 

.-266 

-.212 

8-14 

-.044 

.169 

.330 

.410 

.394 

.288 

.143 

TABLE  XII 

S-ARRAY  at  f = (-l)”'p  for  the  100  Sunspot  Numbers  1770  to  1869 
m m 

m S^(m)  S2(m)  S^Cm)  S^(m) 

-10  -1.8062  3.6869  -6.2717  16.9809 

-9  -1.4952  3.4112  -7.6894  24.5965 

-8  -.7331  3.3526  -9.3130  -45.6024 

-7  -5.8468  3.5271  -11.2804  7.1101 

-6  -2.2570  4.0622  4.5181  -.1710 

-5  -1.6366  4.8009  -30.8892  55.7477 

-4  -.5891  4.7937  2.4144  55.3535* 

-3  -7.1500  4.9124  -39.6295*  3.3841* 

-2  -2.8833  4.6547*  -3.1892*  3.3249 

-1  -2.2403*  2.9516*  -.9613  -.2644 

0 -1.8062*  2.8146  -3.4689  14.3601 

1 -1.5310  2.9172  48.0260  4.7632 

2 -1.1626  2.9305  -5.2351  10.9746 

3 1.4339  3.0070  -5.9414  15.1016 

4 -2.5709  3.3499  -6.8192  -3.6098 

5 -1.7956  3.6912  -8.2632  16.0930 

6 -1.2063  3.8460  -6.0742  -5.0782 

7 2.7461  3.7688  -1.7353  -.2789 

8 -3.0195  3.5825  7.7039 

9 -2.2404  3.7350  -5.4833 

10  -1.9612  3.1968  • 

11  -1.7313  2.8076 

12  -1.4964 

Before  proceeding  to  the  nonstationary  case  we  mention  that  the 
above  technique  can  be  done  graphically  if  so  desired  although 
some  information  would  be  lost.  The  graphical  presentation  of 


Example  1 is  given  below  for  f “ {-l)"p  and  it  easily  identifies 
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IV.  THE  D(n,m)  STATISTIC 


The  identification  procedure  thus  far  suggested  is  one  of 
pattern  recognition.  In  this  section  a statistic  is  suggested  to 
measure  agreement  with  the  proper  pattern  for  a stationary  ARMA 
(P>  <l)  process.  There  is  no  question  that  some  information  will 
be  lost  in  this  statistic.  However,  as  will  be  seen,  it  is  suffi- 
ciently sensitive  to  different  values  of  p and  q that  in  many 
problems  it  will  be  quite  adequate  for  identifying  p and  q.  In 
fact  in  all  of  the  previous  examples  this  statistic,  which  we  refer 
to  as  the  D-Statistic,  easily  identifies  the  "proper"  values  of  p 
and  q.  In  some  later  small  s^unple  exeimples  the  statistic  shows 
appropriate  contending  models  and  suggests  the  obvious  fact  that 
a visual  examination  of  the  S and  R arrays  is  in  order. 

Definition  4. 

For  each  pair  of  integers  (n,m)  such  that  n = 0,  1,  2,  ...,  N 
and  m = 0,  1,  2,  ...,  M,  define 


D (n,m) 


®n+l^^-m-n-l^ 


S (f 
n -m-n) 


,2 


^n^^m-n+1^  ^n+l^^m-n' 


^n-l^^m-n+2^  ®n^^m-n+l^ 


i=l 


'W(^-m-n-i> 


“^n^^-m-n-i+l^ 


Vl<Wi>  V 


*^n  ^^m-n+i+1^ 


where  S , (f  ) =0  and  r (f  ) =1. 

-1  m o m 

Inspection  of  Tables  III  and  IV  as  well  as  the  recursion  rela- 
tions (4)  and  (5)  furnish  the  motivation  for  our  definition  of  D(n,m). 
Although  not  clear  from  Tables  III  and  IV  the  relations  (4)  and  (5) 
show  the  zero  behavior  in  the  R-array  will  tend  to  continue  in  columns 
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p+2,  p+3,  ...  due  to  the  multiplication  by  r in  the  calculation 

n 

of  Also  from  those  relations  it  is  clear  the  "infinity  and 

-Cj^  effect"  will  tend  to  continue  in  columns  p+1,  p+2,  ...  . This 
effect  is  nullified  however  in  D(n,m)  by  the  divisions  in  the  nu- 
merator and  denominator  emd  hence  a large  value  of  D(n,m)  is  sug- 
gestive of  agreement  with  an  ARMA  (n,m)  process.  Thus  the  D- 
estimate  for  p and  q is  defined  by  the  rule,  select  the  model  ARMA 
(P/q)  where 

D(p,q)  = Max  D(n,m) . 

n=l, . . . ,N 
m»0, . . . ,M 

The  case  n = 0 is  excluded  in  the  definition  of  D(p,q)  since  it  is 
not  "normalized".  See  Sections  VII  and  IX. 

Although  these  authors  feel  it  is  always  instructive  to  in- 
spect the  R and  S arrays,  when  D(p,q)  is  as  distinctive  as  it  is 
for  the  examples  just  considered,  the  statistic  will  normally 
suffice  if  stationarity  is  assured. 

However  if  the  process  is  nonstationary  or  seasonal  in  the 
sense  that  columns  in  the  R-arrayhave  runs  of  "near  zeros" interlaced 
among  values  not  near  zero,  the  D-Statistic  would  not  normally  be 
of  value  and  hence  as  the  remaining  part  of  the  paper  will  show  a 


visual  examination  of  the  R and  S arrays  is  usually  better.  We 
now  list  a short  list  of  D(n,m)  values  for  the  previous  examples. 
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In  Section  VII  we  will  investigate  smaller  seunple  cases  andO(n,m) 
will  not  be  nearly  so  distinctive  as  these  examples.  Neverthe- 
less it  will  be  found  to  be  useful  even  in  small  sample  cases. 

A limited  Monte  Carlo  study  comprised  of  approximately  10 
different  processes  and  50  repetitions  of  each  at  a variety  of 
Scunple  sizes  was  performed  by  the  authors  to  compare  the  D- 
Statistic  as  an  automated  method  for  estimating  p and  q with  exist- 
ing methods.  The  methods  compeured  against  were  the  FFE  and  AIC 
methods  suggested  by  Akaike  (1,2)  and  a method  suggested  by 
Anderson  (3).  Of  course  the  FPE  and  Anderson  method  must  al- 
ways select  an  AR(p)  process  but  the  comparisons  were  still  of 
interest.  For  strictly  autoregressive  processes  the  Anderson 
method  always  did  the  best  with  the  D-Statistic  close  behind.  For 
processes  in  general,  i.e.,  ARMA  (p,  q) , the  0-Statistic  performed 
much  better  than  its  only  competitor,  the  AIC  criteria.  In  addi- 
tion to  correctly  identifying  p,q  far  more  often  than  the  AIC 
criteria,  the  computation  time  for  the  D-Statistic  is  a great 
deal  smaller. 

V.  ARMA  (p,  q)  PROCESSES  - THE  NONSTATIONARY  CASE 

In  our  previous  examples  we  have  assumed  a stationary  model. 
For  most  problems  this  is  not  a realistic  assumption  and  some 
attempt  must  be  made  to  transform  the  data  to  stationary  data.  In 
the  case  of  ARMA  (p,  q)  modeling  the  most  common  approach  is  to 
attempt  to  fit  an  ARIMA  (p,  d,  q)  model.  That  is,  observe  the 
sample  autocorrelation,  if  it  fails  to  damp  out  "sufficiently 
fast",  difference  the  data  (possible  several  times)  until  the  re- 
sulting data  produces  an  autocorrelation  which  deimps  out  rapidly. 
Of  course  even  though  the  process  may  be  ARMA  (p,  q)  with  roots 
on  the  unit  circle,  differencing  may  not  produce  a stationary 
process.  In  fact  it  is  clear  that  for  such  models  the  procedure 
will  only  be  beneficial  when  1 is  "almost  a solution"  of  the 
characteristic  equation. 


In  this  section  the  theorems  of  the  previous  sections  for 
model  identification  will  be  generalized  to  the  nonstationary 
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ARMA  (p,  q)  process.  As  a result  it  will  be  shown  how  a transfor- 
mation to  stationarity  may  be  obtained  anytime  the  process  is  non- 
stationary due  to  roots,  real  or  complex,  on  the  unit  circle. 

Recall  that  an  ARMA  (p,  q)  process  is  nonstationary  if  and  only  if 
the  roots  of  the  characteristic  equation 

1 - i^.x  - - ...  - (J  x^  =*  0 (19) 

12  p 

lie  outside  the  unit  circle.  If  the  roots  lie  much  in  the  interior 
of  the  circle  then  the  time  series  is  quite  explosive  and  easily 
recognized  as  nonstationary.  In  this  paper  we  will  not  consider 
such  series  but  will  concentrate  on  nonstationarities  due  to  roots 
on  the  unit  circle.  Thus  in  the  remaining  portion  of  the  paper  our 
use  of  the  word  nonstationary  will  mean  that  the  characteristic 
equation  has  roots  ^n  the  unit  circle.  We  give  the  following  def- 
inition and  theorem  for  future  reference. 

Definition  5. 

An  ARMA  (p,  q)  process  will  be  said  to  be  purely  nonstationary 
if  all  of  the  roots  of  its  characteristic  equation  lie  on  the  unit 
circle. 


Theorem  6. 


Let  r , r , . . . , r be  the  roots  of  the  characteristic  equation 

^ ^ 2 ^ D 

4i(x)  = 1 - iji  X - <(>,x  - ...  - ||)  x^  = 0 of  an  ARMA  (p,  q)  process  and 

12  p 

suppose  the  ' s are  such  that  ~ •••'  P*  Then 

|(fi  I = 1 iff  the  process  is  purely  nonstationary. 

Proof.  Since 

^ n r"^ 

^ 1=1 


we  have 


But  since  |r^|  > 1 for  all  i the  theorem  follows  at  once. 
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In  order  to  interpret  when  the  underlying  process  is  non- 
stationary we  make  the  following  definition. 

Definition  6. 


Denote  by  p (a,,  o., 
m X ^ 


a ) the  autocorrelation  of  an  AKMA 
P 


(p,  q)  process  with  roots  of  the  characteristic  equation  a^,  such 

that  la. I > 1.  Let  r, , r-,  ...  r be  the  roots  of  the  character- 
1 1 ^ P 

istic  equation  of  a given  ARMA  (p,  q)  process  with  j roots  on  che 
unit  circle,  which  for  convenience  we  denote  by  r^^,  •••  tj- 

Then  define 


0* 

m 


and 


where 


lim  ^ 

|a^|,...|a  l-^l 


‘l'“2' 


■j+r 


1,2. . .p. 


2^  1 i'*2  I • • * i^j  I " 

The  condition  in  Definition  6 restricting  the  path  along  which 

the  limit  is  taken  is  necessary  to  guarantee  the  existence  of  the 

limit.  It  is  compatible  with  the  idea  that  in  practice  roots 

"near"  and  equally  close  to  unit  circle  would  be  classified  as  on 

the  unit  circle.  Henceforth  we  let  p = p*  so  that  p is  well  de- 

mm  m 

fined  for  both  stationary  and  nonstationary  processes. 

We  are  now  in  a position  to  show  the  basic  result  which  we 
will  utilize  to  detect  nonstationarity  in  data. 

Theorem  7. 

Let  p = p*  and  assume  Ir. 1 >1.  Suppose  there  exists  an 
mm  ' 1 ' - 


il  > 0 which  is  the  smallest  integer  such  that  S.  (f  ) 

ic  m 


Cg  for 


m = 0 + 1,  +2, 


where  is  a constant  depending  on  whether 


f = p or  f = (-1)  p . 
m m m m 

Proof.  Assume  w.l.o.g.  that  f = p . 
mm 


Then  the  process  is  nonstationary. 


Then  if  (m) 


^ 0 for 


0 + 1,  + 2, 


it  follows  from  Theorem  8 of  [7]  that  p 


, th 


m 


satisfies  an  I order  difference  equation.  More  specifically. 


there  exists  constcuits  , tfi 


1 + 

m 1 m-i 


1'  '^2'  ••• 


such  that 

0,  +1,  +2, 


for  m 


But  (m)  = 3^  0 implies  from  the  proofs  of  Theorems  2 and  3 

that  (p„  = -1  and  from  Theorem  6 p satisfies  the  characteristic 
jC  m 
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equation  of  a purely  nonstationary  process  and  hence  j > 1 in 
Definition  6 and  the  process  is  nonstationary.  If  = 0,  since 
i is  minimal,  (l;k)  0.  However  the  proof  of  Theorem  3 in 

[7]  then  implies  1 - E $.  = 0 which  again  implies  j > 1 so  that 

.-u  . i=l  . 

the  process  is  nonstationary. 

In  order  to  exemplify  Theorem  7,  the  S-Array  at  f = (-l)"'p 

m m 

for  a realization  of  length  200  from  the  AR(2)  process 

(1-.99B)  (1-.5B)X^  = \ - 1.49X^_j^  + .495X^_2  = Z^,  (20) 

where  Z “v  N(0,  1),  followed  by  the  S-Array  at  f = (-l)"*p  , is 
t 'mm 

given  in  Table  XIII.  In  this  case  the  root  of  interest  (call  it 

r^)  is  slightly  outside  the  unit  circle.  If  it  were  exactly  1,  as 

a later  result  will  show,  the  column  would  be  identically  -2. 

Thus  in  that  event  H = 1 in  Theorem  7.  The  process  of  Equation 

(20)  is  of  course  mathematically  stationary.  However,  in  practice 

such  an  equation  is  a nrnperfect  model  and  consequently  for  most 

purposes  we  would  prefer  to  classify  the  model  as  nonstationary 

when  it  has  an  estimated  root  so  near  the  unit  circle.  Thus  if 

such  a model  were  obtained  via  a set  of  data  the  factor  (1-.99B) 

would  be  altered  to  (1-B) . This  will  in  general  be  our  modus 

operand!  for  roots  indicated  near  the  unit  circle,  real  or  complex. 

Moreover  since  our  use  of  the  S-Array  will  therefore  be  to  detect 

an  indication  of  a root  near  the  unit  circle  our  simulated  examples 

will  be  for  "nearly  nonstationary"  processes.  It  should  be  noticed 

that  when  the  sample  autocorrelation  is  utilized  in  Table  XIII 

that  the  constant  behavior  which  is  supposed  to  appear  in  column 

two  if  the  process  is  stationary  has  been  completely  obscured  by 

I the  nonstationary  nature  of  the  data.  The  next  few  theorems  and 

examples  should  make  clear  this  behavior  as  well  as  establish  how 

1 to  deal  with  it.  Before  leaving  this  example  however  let  us  note 

[ that  the  autocorrelation  here  tends  to  a solution  of  a first  order 

I difference  equation  as  r,  ->■  1.  Thus  the  P (r,  ,r.)  satisfies 

I 1 ^l"^^  m 1 2 

I a lower  order  difference  equation  that  the  characteristic 

j equation,  namely  the  one  generated  by  the  nonstationary  factor. 

[ 

I 

I 

I 
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That  is  in  this  case  if  p*  = lim  p (r, ,r^)  we  have 

m m 1 2 

(l-B)p*  « 0,  m = 0,  + 1,  ...  , and  it  is  clear  why  in  Tedale  XIII, 
m 

S^(m)  ^ -2.  Since  in  practice  one  essentially  assumes  the  auto- 
correlation of  a nonstationary  process  is  a limit  such  as  p*, 

m 

Definition  6 and  Theorem  7 are  of  practical  importance. 


m 

Sj^(m) 

TABLE  XIII 

in^ 

f = (-1)  p^ 
m m 

S2(ni) 

S2(m) 

-8 

-2.032 

1.826 

-3.438 

-7 

-2.034 

-5.808 

-5.032 

-6 

-2.034 

4.214 

-52.454 

-5 

-2.034 

2.895 

-21.068 

-4 

-2.035 

-33.755 

11.327 

-3 

-2.035 

5.702 

37.676* 

-2 

-2.031 

16.718* 

-2.115* 

-1 

-2.024* 

2.240* 

-1.233 

0 

-1.976* 

3.053 

-2.542 

1 

-1.969 

1.854 

-4.025 

2 

-1.965 

6.615 

-8.957 

3 

-1.966 

3.803 

6.454 

4 

-1.966 

1.456 

5.604 

5 

-1.967 

-17.359 

-12.754 

6 

-1.966 

-1.906 

-5.647 

7 

-1.968 

4.721 

-2.597 

8 

-1.965 

2.426 

-20.306 

-8 

-2.010 

f = (-l)X 

m m 

6.030 

u 

-7 

-2.010 

6.030 

u 

-6 

-2.010 

6.030 

u 

-5 

-2.010 

6.030 

u 

-4 

-2.009 

6.030 

u 

-3 

-2.008 

6.030 

00* 

-2 

-2.007 

6.030* 

-2.985* 

-1 

-2.003* 

2.985* 

u 

0 

-1.997* 

2.985 

u 

1 

-1.993 

2.985 

u 

2 

-1.992 

2.985 

u 

3 

-1.991 

2.985 

u 

4 

-1.990 

2.985 

u 

5 

-1.990 

2.985 

u 

6 

-1.990 

2.985 

u 

7 

-1.990 

2.985 

u 

8 

-1.990 

2.985 

u 
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Theorem  8. 

Suppose  an  ARMA  (p,  q)  process  has  k roots  of  the  character- 
istic equation  which  approach  the  unit  circle  as  in  Definition  6, 

N of  which  are  distinct.  Then  for  m = 0,  +1,  p*  satisfies 

m 

a linear  homogeneous  difference  equation  of  order  less  than  or 
equal  to  N,  whose  characteristic  equation  has  all  of  its  roots 
distinct  and  on  the  unit  circle. 

Proof.  The  proof  of  Theorem  8 is  included  in  Appendix  3 for  an 
ARMA  (3,  1)  process  only.  It  was  first  shown  by  Kelley  in  [8] 
where  the  cases  ARMA  (1,  0),  (1,  1),  (1,  2),  (2,  0),  (2,  1),  (2,  2) 
and  (3,  0)  are  also  given.  In  addition  the  result  is  demonstrated 
for  an  ARMA  (4,  1)  in  Example  6.  An  analytic  proof  of  higher  order 
cases  seems  to  be  intractable/  at  least  by  Kelley's  approach,  due 
to  the  lengthy  algebra  involved.  However  since  Theorem  8 and 
Theorem  9 which  follow  will  only  be  utilized  to  suggest  an  appro- 
priate transformation  to  stationarity , we  shall  be  satisfied  with 
•cne  validity  of  the  result  as  given.  That  is,  once  the  transfor- 
mation is  constructed,  it  can  be  observed  whether  or  not  it  had 
the  desired  effect.  As  our  final  theorem  we  make  us  of  Theorem  8 
to  obtain  the  converse  of  Theorem  7. 

Theorem  9. 

Suppose  Theorem  8 holds  then  the  condition  that  (m)  is 
constant  for  some  J,  < p is  both  necessary  and  sufficient  for  non- 
Etationarity. 

Proof.  The  sufficiency  follows  of  course  from  Theorem  7.  Now  let 

^m  ~ *^m  suppose  the  process  is  nonstationary  with  autocorrelation 

p = p*.  Then  let  N be  the  number  of  distinct  roots  approaching  the 
m m 

circle  in  the  manner  described  by  Definition  6.  By  Theorem  8,  p 
satisfies  the  characteristic  equation  of  a purely  nonstationary  AR 
process  of  order  i < N,  i.e. 

(1  - - ...  - \().B**’)p  = 0 (21) 

i z ic.  m 


and 


But  since  satisfies  an  equation  such  as  (21),  1 - 0, 

the  proof  of  Theorem  1 (see  [6,7])  implies  Theorem  2 and  ^ Theorem  3 
hold  even  though  the  process  is  nonstationary. 


Thus  = -1  implies 


or  ()c)  constant.  On  the  other 


hand  if  ip.  = 1,  since  all  the  roots  are  on  the  unit  circle,  it  is 

easily  shown  that  1 - = 0.  However  in  this  case  the  proof 

given  of  Theorem  1 in  [ 6 ] can  be  modified  to  show  (ra)  = 0 so  that 

again  we  have  S.  (m)  constant.  The  proof  for  f = (-l)"'p  is 
io  in  m 

entirely  similar  and  hence  is  not  argued.  Before  proceeding  we  list 
three  useful  corollaries  to  the  above  results  whose  proofs  follow 
from  the  proofs  of  those  results. 

Corollary  2.  If  S, (m)  = 0 when  f = p or  if  S, (m)  = -2  when 
1 m m 1 


istic  equation  has  at  least  one  root  of  +1  on  the  unit  circle. 

If  S, (m)  = -2  when  f = p or  S, (m)  = 0 when 
1 m m 1 


B 

II 

(-dX' 

m 

istic 

equation 

Corollary  3. 

f = 
m 

(-1)  P , 
m 

istic 

equation 

Corollary  4. 

(-1)  P 


then  the 


2 mm  m m 

process  is  nonstationary  and  the  characteristic  equation  has  at 

least  one  root  of  +1  and  at  least  one  root  of  -1.  A number  of  re- 
sults such  as  Corollaries  2,  3 and  4 could  be  listed.  They  of  course 
demonstrate  the  manner  in  which  evidence  from  the  S-Array  can  be 
useful  in  determining  the  appropriate  transformation  to  stationarity . 
We  now  list  a final  corollary  of  this  nature.  Others  will  become 
apparent  in  the  next  section. 
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Corollary  5.  If  (m)  = C 0,  then  the  process  is  nonstationary 
and  the  characteristic  equation  has  at  least  two  complex  roots  on 
the  unit  circle.  Tables  XIV  and  XV  illustrate  the  above  results  as 
a root  approaches  1^  for  an  ARMA  (1,  1)  and  an  ARMA  (3,  1).  In  both 
cases  the  convergence  of  (k)  to  -2  is  clear  even  though  the>  sample 
autocorrelations  were  utilized  in  the  calculations.  Example  6 which 
follows  demonstrates  Corollary  5.  Theorem  7 is  demonstrated  graphi- 
cally in  Figures  1,  2 and  3.  It  should  be  noted  that  although  the 
autocorrelations  are  dramatically  different,  the  S^(m)  columns  are 
essentially  the  same.  That  is,  they  are  roughly  constamt  with  pos- 
sibly some  slight  deviation  about  the  "center  line",  k = -£  + ^.  The 
reason  for  the  possible  deviation  can  be  seen  from  Tables  XIV  and  XV 
where  it  is  clear  that  the  elements  S^(-l)  and  Sj^(O),  due  to  the  mov- 
ing average,  behave  somewhat  differently  than  the  remaining  elements 
although  that  difference,  as  it  must,  is  dimensioning  as  -*•  1. 

Let  us  now  consider  some  additional  examples  which  make  use  of 
all  of  our  results  to  this  point. 

VI . APPLICATIONS 

The  examples  which  follow  include  some  from  generated  data  and 
some  from  real  data.  In  the  case  of  generated  data,  the  data  (as  in 
the  previous  examples)  was  generated  on  the  CDC  Cyber  72  via  the 
ARMA  (p,  q)  equation  with  Gaussian  white  noise  input  unless  stated 
to  the  contrary . 

Example  5.  (Series  C) 

The  source  of  data  for  this  example  are  226  temperature  readings 
from  a chemical  process  taken  every  minute.  The  data  is  given  in 
[ 5 ] where  it  is  referred  to  as  Series  C.  The  data  is  plotted 
below  in  Figure  4. 
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Tables  XVI  and  XVII  show  the  S-array  for  the  data.  Table  XVI 
clearly  identifies  the  process  as  essentially  nonstationary  due  to 
a unit  root.  In  addition  that  array  seems  to  clearly  indicate  the 
process  is  aji  AKMA  (2,  1).  However,  a word  of  caution  is  in  order. 
That  is,  as  we  have  previously  seen  and  explained,  the  appcirent 
nonstationarity  will  often  completely  dominate  the  column  behavior. 
In  Tedsle  XVI  it  is  not  clear  that  this  is  the  case  and  the  ARMA 
(2,  1)  shows  through  remarkedly  well.  Nevertheless  the  non- 
stationarity should  be  removed  before  the  final  identification. 
Table  XVII  shows  the  S-array  for  the  first  differences  of  Series  C 
The  differences  are  by  that  table  a stationary  AR  (1)  with 


1.805  > 
2.241 


.81 


Thus  the  model  obtained  is 

(1-B) (1-.81B)  . (22) 

It  should  now  be  clear  that  the  apparent  moving  average  term  in- 
dicated in  Table  XVI  is  indeed  induced  by  the  nonstationarity 
which  is  vividly  displayed  in  Column  1 of  that  tadsle.  As  previously 

explained  this  is  to  be  expected  since  p*  corresponding  to  (22) 

m 

actually  satisfies  a first  order  DE.  The  remarkable  thing  here  is 
that  the  2nd  order  process  shows  up  at  all  in  Table  XVI.  Before 
leaving  this  example  let  us  make  a point  which  will  be  useful  in 
the  next  example.  That  is,  notice  that  we  have  actually  obtained 
estimates  of  the  coefficients  one  factor  at  a tine.  From  S^(m)  of 
Table  XVI  we  obtained  the  factor  (1-B)  and  from  Table  XVII  we 
obtained  the  estimate  The  resulting  model  is  the  same  as  the 
one  given  for  this  data  in  Box  and  Jenkins  [5] . The  next 
example  shows  how  this  observation  can  be  extended. 

Example  6. 

We  now  consider  an  example  of  an  ARMA  (4,  1)  process  with  two 
complex  roots  nearly  on  the  unit  circle.  The  data  consists  of  300 
points  generated  from  the  process 

(1-1.2B+.45B^) (1-1.7B+.99B^)  = (1-.5S) 


(23) 
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TABLE  XVI 

A Portion  of  the  S-Array  at 
fm  = 

for  the  Series  C Data 


TABLE  XVII 

A Portion  of  the  S-Array  at 
f = C-D^P 


m m 

for  the  First  Difference 
of  the  Series  C Data 


m 

§1 

^2 

§3 

m 

§1 

^2 

-7 

-2.080 

4.755 

-2.464 

-7 

-2.217 

5.250 

-6 

-2.073 

4.314 

-5.288 

-6 

-2.192 

-.212 

-5 

-2.064 

4.471 

5.590 

-5 

-2.163 

7.000 

-4 

-2.056 

4.497 

567.959 

-4 

-2.190 

1.911 

-3 

-2.046 

4.449 

13.380* 

-3 

-2.240 

-.969 

-2 

-2.035 

9.576* 

-2.100* 

-2 

-2.234 

-170.077* 

-1 

-2.022* 

2.492* 

-3.544 

-1 

-2.241* 

1.786* 

0 

-1.977* 

3.604 

-1.989 

0 

-1.805* 

.546 

1 

-1.965 

a.572 

-23.995 

1 

-1.810 

-10.704 

2 

-1.955 

3.589 

4.549 

2 

-1.806 

2.706 

3 

-1.946 

3.704 

2.909 

3 

-1.840 

.164 

4 

-1.939 

3.413 

-1.719 

4 

-1.859 

3.127 

5 

-1.931 

4.060 

-4.178 

5 

-1.838 

11.435 

6 

-1.925 

3.660 

31.451 

6 

-1.821 

1.314 

7 

-1.918 

3.604 

-2.994 

7 

-1.712 

12.830 

A graph  of  the  data  is  given  in  Figure  5 while  a graph  of  the  auto- 
correlations is  given  in  Figure  6.  In  Table  XVIII  the  S-Array  for 
this  data  is  shown  out  to  Column  5.  The  fact  that  (m)  is  nearly 
a nonzero  constant  immediately  implies  the  process  has  two  complex 
roots  very  near  the  unit  circle  and  that  the  sample  autocorrelation 
behaves  as  the  autocorrelation  of  a purely  nonstationary  AR  (2)  pro- 
cess with  ^ Thus  we  assume  the  process  to  have  the  factor 

1 - 0j^B  + But  since  the  seunple  acf  behaves  as  that  of  an  AR  (2) 

we  have  from  our  previous  results  that 

1 + ^ 

= C - 2 

Now  if  we  estimate  C by  averaging  several  of  the  elements  in  Column  2 
surrounding  the  starred  values,  say  S2  (-4)  to  S2(l),  we  obtain  the 
estimate  C = 3.71  and  hence 

= 1.71 
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TABLE  XVIII 

jn** 

A Portion  of  the  S-Array  at  f = C-1)  P for  a Realization  of  Length 

mm 

300  of  the  ARMA(4,1)  Process  in  Example  6 


k 

§1 

§2 

§3 

^4 

^5 

-8 

-3.177 

3.736 

-5.117 

4.003 

-27.408 

-7 

-2.257 

3.734 

-12.166 

28.431 

220.437 

-6 

-1.917 

3.733 

-7.187 

26.320 

91.826 

-5 

-1.619 

3.732 

2.208 

24.901 

81.516* 

-4 

-1.089 

3.732 

-21.221 

61.184* 

-5.304* 

-3 

8.634 

3.731 

-12.892* 

5.673* 

-5.572 

-2 

-2.825 

3.724* 

-5.191* 

9.285 

-6.364 

-1 

-2.167* 

3.701* 

-4.477 

7.083 

-179.415 

0 

-1.856* 

3.691 

-1.372 

7.124 

-5.327 

1 

-1.547 

3.690 

-7.677 

-6.694 

-6.548 

2 

-.896 

3.690 

-5.323 

3.634 

-8.502 

3 

-12.228 

3.690 

-13.650 

9.761 

2.505 

4 

-2.613 

3.689 

-7.846 

1.059 

-13.610 

5 

-2.089 

3.688 

-7.124 

-2.673 

-49.582 

6 

-1.794 

3.686 

-7.723 

39.953 

4.982 

7 

-1.459 

3.685 

-6.926 

7.759 

-11.612 

The  operator  1-1. 

71B  + B^ 

thus 

constructed 

is  a nonstationary 

operator 

and,  unless  it  is  a 

repeated  factor. 

Table  XVIII 

suggests 

that  the 

data  transformed  by 

this 

operator  will  be  stationary.  Thus 

let  w = 
u 

(1-1.71B+B 

)X^.  Table  XIX  gives  the 

R and  S-Arrays  for  w^. 

Clearly 

from  (m) , 

p = 2,  q 

= 1. 

The  fact  that  q = 1 is 

also  clear 

in  (m) . The  model  thus  obtained  is 

= (l-eB)Z^,  (24) 


where  ({i  (B)  is  a stationary  operator.  Now  from  the  IMSL  subroutine 
FTMAXL  one  obtains  from  the  w^  series  (shown  in  Figure  7 on  an  ex- 
panded scale) 

ij>(B)  = 1 - 1.17B  + .45B^ 

0 = .52 

and  hence  we  finally  obtain  the  model 

(1-1.17B+.45B^) (1-1.71B+B^) (X^-X)  = (1-.52B)Z^,  (25) 

which  is  very  close  to  (23)  . In  fact  if  we  agree  the  data  should  be 
modeled  as  nonstationary,  this  model  is  about  as  good  as  one  could 
expect  to  get. 


I 
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TABLE  XIX 


f 

[ 

i 

i 


! 


Portions  of  the  R and  S 
2 

= (1-1.71B+B  )X^, 


R-Array 

k 

5i 

^2 

-8 

-.195 

.001 

-7 

.186 

.023 

-6 

-.181 

-.016 

-5 

.129 

.047 

-4 

-.065 

-.683 

-3 

-.037 

-.156 

-2 

.264 

.112 

-1 

-.601 

-.398* 

0 

1.000* 

.060 

1 

-.601 

-.054 

2 

.264 

.061 

3 

-.037 

.335 

4 

-.065 

-.025 

5 

.129 

.027 

6 

-.181 

.001 

7 

.186 

.019 

8 

-.195 

-.037 

Arrays  at 

f = 
m 

C-D^P 

m 

where  X^ 

is  the 

Series 

?3 

k 

§1 

-.026 

-8 

-1.952 

-.016 

-7 

-1.970 

.005 

-6 

-1.717 

.036 

-5 

-1.503 

.013 

-4 

-.419 

.129 

-3 

-7.971 

.337* 

-2 

-3.272 

.029 

-1 

-2.661’ 

.005 

0 

-1.601’ 

.019 

1 

-1.440 

.002 

2 

-1.143 

-.013 

3 

.722 

-.017 

4 

-2.987 

-.020 

5 

-2.394 

-.024 

6 

-2.030 

.015 

7 

-2.049 

-.009 

8 

-1.849 

for  the  Series 
of  Figure  5 


S-Array 

S., 

S, 

-2 

-3 

-24.780 

-1.064 

2.909 

-7.655 

5.878 

41.254 

6.501 

-3.918 

6.152 

48.943 

5.627 

19.511* 

12.102* 

-1.686* 

1.846* 

-2.190 

2.732 

5.600 

2.442 

-2.778 

3.193 

-21.343 

3.216 

1.094 

4.953 

.390 

1.898 

-3.788 

-19.886 

-8.893 

5.495 

-5.360 

3.294 

-.412 

Example  7. 

Table  XX  gives  the  S-Array  at  f = p to  4 Columns  from  a real- 

m m 

ization  of  length  300  froa  an  AR  process.  Since  S^im)  is  essentially 
constant  the  process  is  nonstationary  with  at  least  three  roots  on 
the  unit  circle.  Since  S^Cm)  ^ 0 , +1  is  not  a root.  Thus  we  have 
at  least  one  root  of  -1  and  two  complex  roots  on  the  unit  circle. 

Hence  if  we  first  let  evaluate  the  S-Array  for 

w^,  the  complex  operator  for  stationarity  is  easily  found  as  in  the 
previous  example  and  the  final  factor  is  obtained  from  the  resulting 
stationary  series.  The  actual  process  for  this  realization  is 

(1-G  B)  (1-G  B)  (1-G  B)  (1-G  B)X  = Z , (26)  j 

where  i z J 4 t t ^ 

Gj^  = -.98,  G^  = .6,  Gj  = .7  + .7i,  = . 7 - .7i 

and  hence  the  validity  of  our  analysis  is  clear. 
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TABLE  XX 

A Portion  of  the  S-Array  at  f = p for  a Realization  of  Length  300 

m m 

4 

from  the  AR(4)  Process 


= -.98, 

= .6,  G3  = . 

.7+.7j,  G^ 

li 

1 

k 

s. 

S- 

S, 

s. 

-1 

-2 

-3 

-4 

-10 

8.8240 

.6895 

-1.2095 

.5710 

-9 

.3983 

.5937 

-1.1984 

.2714 

-8 

-.3445 

.4755 

-1.2022 

2.3810 

-7 

-1.0536 

.7379 

-1.1975 

-3.9603 

-6 

21.1401 

.4903 

-1.1944 

-.7487 

-5 

.3250 

.5849 

-1.2075 

-3.5518 

-4 

-.2808 

.6935 

-1.1784 

-1.5304* 

-3 

-1.0689 

.4218 

-1.2420* 

.6201* 

-2 

14.2091 

.7097* 

-1.0425* 

.8323 

-1 

.4601* 

.5667* 

-1.1233 

.4306 

0 

-.3151* 

.4468 

-1.1084 

.8622 

1 

-.9342 

.7213 

-1.1252 

2.2520 

2 

-15.4989 

.4718 

-1.1223 

-.3280 

3 

.3904 

.5741 

-1.1193 

-1.0177 

4 

-.2453 

.6628 

-1.1249 

-2.1680 

5 

-.9548 

.4323 

-1.1182 

.6945 

6 

-19.6517 

.6786 

-1.1486 

.8472 

Example  8. 

Figure 

8 below  shows 

a graph  of 

the  sample 

acf  of  a realiza- 

tion  of  length  300  from  an  ARtlA  process. 


Figure  8 


R(K) 


Table 

XXI  shows 

the  S-array 

at  f = (-1) 
m 

TABLE  XXI 

^p (m)  for  the 

data. 

m 

S, 

S, 

S, 

1 

2 

3 

4 

5 

-6 

-1.901 

.649 

-3.211 

12.855 

-5.575 

-5 

-2.389 

4.077 

-3.010 

7.945 

60.297* 

-4 

-2.068 

-.627 

-3.116 

8.529* 

-4.163* 

-3 

-1.815 

2.492 

-3.795* 

4.472* 

8.970 

-2 

-2.181 

-19.321* 

-2.933* 

4.752 

-1.963 

-1 

-2.234* 

1.654* 

-2.481 

3.812 

-6.840 

0 

-1.809* 

17.868 

-2.888 

4.954 

-8.589 

1 

-1.846 

.497 

-2.919 

5.765 

-9.296 

2 

-2.227 

3.489 

-2.783 

6.196 

50.935 

3 

-1.936 

-.787 

-2.814 

6.389 

-5.375 

4 

-1.719 

2.215 

-2.951 

4.856 

-25.993 

5 

-2.109 

-43.331 

-2.794 

4.687 

-7.443 

An  inspection  of  the  table  immediately  suggests  from  Column  1 

that  the  data  needs  differenced.  Moreover  from  Columns  3 and  4 one 

would  assume  (by  Column  3)  that  the  process  also  has  two  complex 

roots  on  the  unit  circle  and  is  a fourth  order  process.  However, 

before  making  decisions  regarding  Colums  3 and  4 one  should  first 

eliminate  the  simplest  indicated  nonstationarity . That  is,  the  data 

should  be  differenced.  Upon  differencing  the  data  one  obtains  Table 

XXII  taking  f = p which  in  this  case  yielded  a more  obvious  pattern 
m « m 

than  f = (-l)''’p  . 
m m 

TABLE  XXII 

S-ARRAY  AT  f = p FOR  (l-B)X^  IN  EXAMPLE  8 
mm  t 


m 

^1 

"2 

-5 

-1.573 

2.535 

1.787 

-4 

.331 

2.522 

2.088 

-3 

-2.275 

2.727 

1.085* 

-2 

-.753 

2.578* 

-.730* 

-1 

-6.103* 

2.230* 

-1.103 

0 

-1.195* 

2.464 

-1.006 

1 

3.050 

2.543 

-1.297 

2 

-1.783 

2.427 

-.927 

3 

-.248 

2.433 

.047 

4 

-2.744 

2.465 

.632 
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Inspection  of  Table  XXII,  then  suggests  the  process  has  two 

complex  roots  on  the  unit  circle.  Proceeding  now  as  in  Example  6 

2 

gives  at  once  the  operator  1 + .51B  + B as  a nonstationary  factor. 


Applying 

this  operator  to 

the 

differenced  data  then  yields  the  R and 

S Arrays 

of  Table  XXIII. 

TABLE  XXIII 

S AND  R-ARRAYS  FOR 

(1+ 

.51B+B^) (1 

-B)X  IN  EXAMPLE  8 

(m) 

S2(m) 

m 

(m) 

Rj  (m) 

-2.494 

.632 

-5 

-.102 

.013 

-2.718 

-1.518 

-4 

.152 

.010 

-2.609 

-11.483 

-3 

-.262 

.016 

-2.712 

9.999* 

-2 

.421 

.087 

-2.385* 

2.080* 

-1 

-.721 

-.278* 

-1.721* 

1.311 

0 

1.000* 

.057 

-1.584 

.581 

1 

.721 

.009 

-1.621 

-.505 

2 

.421 

.006 

-1.581 

. 2.891 

3 

-.262 

.008 

-1.669 

5.530 

4 

.152 

-.006 

-1.770 

-.332 

5 

-.102 

.013 

Finally 

by  the  D-Statistic 

or 

visually  analyzing  Table 

XXIII  one  can 

see  the 

resulting  data  is 

ARMA 

(1,1)  and 

the  final  model  is 

(1-.56B) (1+. 

51B+B^) (l-B)X^ 

= (1+.4B)Z^. 

The  coefficient  .56  is  obtained  from  Table  XXIII  from 

the  ratio 

1.584 

-2) 

2.712 

• DO 

and  as  mentioned  before  is 

the 

conditional  Yule-Wallcer 

estimate  (i.e 

given  p 

and  q) . The  estimate 

. 4 on  the 

right  side  of 

the  equation 

was  obtained  by  standard  methods.  Since  the  data  of  this  example  was 
generated  from  the  model 

(1-.6B) (1+.5B+.98B^) (1-.98B)X^  = (1+.4B)Z^ 
the  above  analysis  appears  quite  good. 
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Example  9. 

Before  returning  to  real  data  let  us  consider  one  more  set  of 
simulated  data.  Again  we  consider  an  AR(4)  process.  Figure  9 shows 
300  generated  points  from  the  process 


(1-G^B) (l-G^B) (l-G^B) (l-G^B)X^  = Z^,  (27) 

where 

G^  = -.99,  “ ■.5+.4i,  Gg  » .5-.4i. 

A 

Table  XXIV  shows  the  S-Array  at  f * p cmd  f =«  (-1)  p . 

m m m m 
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TABLE  XXIV 


^m  ~ ^m 

S-Arrays 

m 

S]^(m) 

S2(m) 

S3(m) 

m 

Sj^(m) 

S2(m) 

,266 

-.073 

.015 

-7 

-2.266 

-.045 

.033 

-.149 

-.077 

.110 

-6 

-1.850 

.004 

.039 

,259 

-.082 

-.726 

-5 

-2,259 

-.037 

.419 

-.136 

-.080 

-.165 

-4 

-1.863 

-.052 

.055 

.249 

-.064 

-.074* 

-3 

-2.249 

.094 

.188* 

-.131 

-.033* 

-.045* 

-2 

-1.868 

-.346* 

-.114* 

.214* 

.028* 

-.054 

-1 

-2.214* 

.291* 

-.018 

-.176* 

.065 

-.069 

0 

-1.S23* 

-.095 

-.040 

.151 

.075 

-.286 

1 

-2.151 

.049 

.101 

-.199 

.077 

.018 

2 

-1.800 

.035 

-.037 

.157 

.074 

-.057 

3 

-2.157 

-.004 

-.042 

-.206 

.069 

-.393 

4 

-1.793 

.042 

.166 

.175 

.068 

.074 

5 

-2.175 

.035 

-.074 

-.210 

.074 

.016 

6 

-1.789 

-.001 

-.075 

.183 

.079 

-.013 

7 

-2.183 

.076 

-.046 

Since  S.. 

2 

(m)  ^ 0 in 

the  interval 

surrounding  the 

starred  values  at 

both  f = p and  f = (-l)"'p  the  process  has  a root  of  both  +1  and 
m m m m 

2 

hence  we  operate  on  the  data  by  the  operator  1 - B . Table  XXV 


shows  the  result  of  that  operation.  From  there  it  is  clear  that 
2 

(1-B  )X  = w is  an  AR(2) . Thus  the  model 

^ 2 2 
(1-i|)^B-,|/2B  ) (1-B 

is  at  once  obtained  where  the  first  factor  is  a stationary  operator. 


Again  the  Yule-Walker  estimate  of  can  be  obtained  from 
Column  2 of  the  S-array.  That  is 


1 - 2,439 

'*'2  “ “ 5.690 


.429  . 


(28) 


Moreover  it  can  be  shown,  as  might  be  guessed  from  Theorem  4, 
that  the  Yule-Walker  estimate,  for  1);^^  can  also  be  obtained 

from  $2 (m) . That  is 

1 + ° 2.439 

il  = 


or 


1.01 


(29) 
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TABLE  XXV 


S-Array  for  w^  = 

(l-B^)X^  at 

m m 

in  Example  9 

i 

i 

S^(m) 

S2(m) 

S2(m) 

F 

i 

-12.445 

3.802 

-221.284 

-6 

-2.438 

3.875 

-1.360 

-5 

-1.675 

7.011 

-10.735 

S -4 

-.655 

6.539 

9.207 

-3 

-11.885 

7.166 

-27.116* 

[ -2 

-3.471 

5.690* 

-2.680* 

! -1 

-2.413* 

2.439* 

-1.308 

i 0 

-1.707* 

2.117 

-7.002 

1 

-1.404 

2.327 

.719 

' 2 

-1.091 

2.737 

-3.084 

3 

1.902 

2.986 

-12.297 

; 4 

-2.481 

3.030 

2.067 

i 6 

-1.695 

3.202 

.729 

-1.087 

3.048 

-8.949 

7 

7.595 

2.918 

-3.245 

Thus  the  estimated  model  is 

(1-1.01B+.429B^) (l-B^)X^  = 
whereas  the  true  model  is 

(1-B+.41B^) (1-.97B^)X^  = . 

Of  course  the  estimates  (28)  and  (29)  can  be  obtained  by  more 
sophisticated  methods  but  in  this  example  its  not  necessary. 
Example  10.  (Series  A) 

Consider  now  the  197  data  points  given  in  Box  and  Jenkins  and 
referred  to  as  Series  A.  The  data  consists  of  concentration  read- 
ings taken  every  two  hours  from  a chemical  process.  A graph  of 
the  data  (connected  by  straight  lines)  is  given  in  Figure  10.  The 
data  are  shown  in  tabular  form  in  [ 5 ] where  it  was  modeled  as 
either  the  stationary  ARMA  (1,  1)  model 

X^.  - •92X^_^  = 1.45  + a^  - •58a^_j^  , (30) 

or  the  nonstationary  ARIMA  (0,  1,  1)  model 
(l-B)X^  “ \ ■ •’^°^t-l  ‘ 

In  Tables  XXVI  and  XXVII  the  S-array  at  f = p and  f = (-1)  p 

mm.  m m 

are  given.  At  f = (-1)  p there  is  clear  support  in  Table  XXVII 
m m 


for  the  models  (30)  and  (31).  However  we  should  remember  that 
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(m)  can  take  on  the  characteristics  shown  there  strictly  because 
of  the  near  or  actual  nonstationarity . Therefore  a close  inspec- 
tion of  the  S-arrays  is  in  order.  Table  XXVI,  although  supportive 
of  (30)  and  (31)  is  not  nearly  so  distinct  as  Table  XXVII  in  t.his 
regard.  Note  that  there  is  a very  distinctive  behavior  in  Table 
XXVI  in  Columns  7 and  8.  That  is,  from  those  columns  the  process 
appears  to  be  either  an  ARMA  (7,  1)  or  ARMA  (7,  0),  and  another 
inspection  of  Table  XXVII  supports  this.  Now  as  stressed  before, 
the  decision  to  include  the  moving  average  term  when  modeling  the 
process  as  a nonstationary  model  should  be  made  after  the  non- 
stationarity is  removed.  Moreover  the  only  thing  that  Table  XXVI 
clearly  indicates  is  that  the  data  should  be  differenced.  Table 

XXVIII  shows  the  S-array  of  the  differenced  data  at  f = p cuid 

m m 

the  process  is  clearly  ARMA  (6,  0) . An  R-array  Table  of  the  data 
is  equally  distinct.  A preliminary  estimate  of  the  coefficients 
in  the  models  (corresponding  to  (30)  and  (31)) 

(X^-X)  (32) 

and 

(1_^  B-ij;_B^-’J;  B^-0^B%.B-i((-B®)  (l-B)X  = a^  (33) 

gives  the  Yule-Walker  estimates 

'•’l  " '*’2  = -2  = -02  <j)^  = .01  = -.01  0^  = .06  0^  = .16 

. 0,  =-.61  0.  = -.40  0_  = -.36  0,  = -.32  0^  = -.31  0,  = -.21  . 

1 J 4 b o 

These  sets  of  estimates  are  compatible  and  suggest  the  parsimonious 
models, 

(1-0^B-02B^-0.^b'^)  (X^-X)  = a^-  9 a^_j_  (34) 

or  2 3 4 5 6 

(1-0^B-02B  -0^3  -0^6  -0^6  -02B  ) (l-B)X^  = a^  . (35) 

Note  that  like  (32)  and  (33),  (34)  and  (35)  are  compatible,  i.e., 

0-  = 0-  = 0.  = 0c  = 0,  implies  0,  = 0.  = 0^  = 0,  = 0.  In  addition 
it  should  be  noted  that  the  nonstationary  ARIMA  (6,  1,  0)  in  (35)  is 
parsimonious  in  that  there  are  only  two  parameters  to  estimate.  Any 
one  of  the  models  (32),  (33),  (34)  or  (35)  yield  a slightly  smaller 
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estimate  of  the  white  noise  variance  than  (30)  or  (31)  . This  is  not 
particularly  significant  inview  of  the  models  even  though  (34)  and 
(35)  are  reasonably  parsimonious. 

There  are  however  a number  of  observations  which  can  be  made 
concerning  Series  A which  suggest  that  the  models  (34)  or  (35)  are 
more  appropriate  models  for  the  data.  In  particular  a spectral 
analysis  of  the  data  by  four  different  methods,  standard  window 

methods,  the  AR  Spectral  Estimator,  the  Maximum  Entropy  Method,  and 
the  G-Spectral  Estimator  all  show  a maximum  pea)c  in  the  spectrum  at 
the  frequency  w = 0 as  well  as  a local  peak  near  w = .12,  while  all 
but  the  windows  show  another  small  peak  near  the  w = .26.  This  of 
course  suggest  the  autocorrelation  is  not  that  of  a process  such  as 
(30)  and  (31).  The  logs  of  the  Parzen  window,  G-Spectral  and  AR- 
Spectral  estimators  are  shown  in  Figure  (11) . Where  the 
log  spectral  density  associated  with  the  model  defined  by  (30)  is 
also  given.  In  that  figure  the  log  AR-Spectral  estimator  can  essen- 
tially be  considered  the  log  of  the  spectral  density  associated 
with  the  model  of  (34)  with  0 =0.  That  is,  in  obtaining  the  AR- 

Spectral  estimator  the  method  of  Akaike  was  utilized  to  obtain  the 
proper  order  underlying  autoregressive  model.  The  model  selected 
was  an  AR  (7,  0)  in  agreement  with  the  order  of  autoregressive  in  (34), 
and  the  estimated  coefficients  were  essentially  the  ip's  already 
given.  Thus  it  is  clear  that  the  estimated  spectral  density  obtained 
in  a variety  of  ways  suggests  the  model  (34)  and  (35)  as  more  appro- 
priate than  (30)  or  (31) . In  addition  that  analysis  suggests  the 
autocorrelation  very  nearly  has  a damped  periodic  component  whose 
period  is  in  the  neighborhood  of  7 or  8 lags.  It  is  quite  possible 
that  such  observations  as  this  may  be  extremely  important  to  the 
experimenter  and  they  are  completely  lost  in  both  of  the  models  (30) 
and  (31).  Of  course,  there  is  no  absolute  answer  here  as  to  whether 
the  models  of  (34)  and  (35)  are  better  than  those  of  (30)  and  (31) 
since  the  true  model  is  not  known.  To  shed  some  additional  light 
on  the  methods  of  this  paper,  the  reasons  why  a model  such  as  (34) 


55 


or  (35)  might  be  identified  as  (30)  or  (31)  or  vice-versa  were 
Investigated  by  a short  Monte  Carlo  study.  The  result  of  simulat- 
ing 5 realizations  of  (30)  of  length  197  was  that  the  AFMA  (1,1) 
model  was  always  selected  via  the  R and  S array  technique.  On  the 
other  hand  when  realizations  of  (34)  of  the  same  length  were 
simulated  the  ARMA  (7,0)  or  ARMA  (7,1)  was  always  chosen. 
Typical  S arrays  for  the  simulated  ARMA  (1,1)  and  ARMA  (7,1)  are 
shown  in  Tables  XXIX  emd  XXX.  From  these  tables  the  ARMA  (1,1) 
model  of  (30)  is  clearly  distinguishable  from  the  model  of  (34) , 
and  hence  these  authors  feel  that  (34)  or  (35)  is  more  appropriate 
for  series  A and  nearly  as  parsimonious  as  the  (1,1)  model. 


Fig.  11  Spectral  Estimators  for  Series  A 
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VII.  SMALL  SAMPLE  SIZES 

In  the  previous  excunples  the  number  of  observations  was  always 
100  or  more.  This  was  done  so  that  patterns  would  be  very  distinct 
to  the  reader  attempting  to  understand  the  method.  However,  once 
the  method  is  understood,  processes  which  are  reasonably  modeled 
by  ARMA  models  can  be  identified  on  much  smaller  sample  sizes. 
Although  it  is  realized  that  the  number  of  samples  required  to 
accurately  identify  p and  q would  vary  dramatically  with  the  coeffi- 
cients, in  order  to  obtain  some  idea  of  how  small  the  sample  size 
might  be  ajid  still  yield  clear  identification  of  p and  q,  the 
authors  generated  3 realization  of  each  of  the  processes  in  the 
previous  examples  at  sample  sizes  25,  50,  75,  100,  150  and  200. 

In  all  but  one  case,  Example  6 which  is  discussed  below,  p was 
clearly  correctly  identifiable  in  sample  sizes  as  small  as  25  in  at 
least  two  of  the  three  realizations.  In  most  cases  q was  either 
correctly  identified  or  identified  as  q + 1 or  q + 2 in  samples  of 
size  25.  The  patterns  in  the  R & S array  were  not  markedly  better 
for  samples  of  size  50  than  for  samples  of  size  25.  However,  by 
75  or  100  samples  the  patterns  were  usually  very  good  and  with  a 
few  exceptions  all  the  realizations  were  correctly  identified, 
even  by  the  D-Statistic.  Although  the  D-Statistic  is  not  recom- 
mended alone  for  small  samples,  these  authors  found  it  very  useful 
even  then.  That  is,  in  each  realization  D(n,m)  was  calculated 
with  Max  n = 4 and  Max  m = 2.  By  investigating  the  R and  S arrays 
at  all  large  values  of  D(n,m)  the  most  predominate  patterns  were 
considered. 

In  order  to  substantiate  the  claims  made  in  the  above,  three 
examples  are  now  given. 

Example  11. 

Let  us  consider  once  more  the  process 

- 1.32X.  , + .68X.  _ = Z.  - .8Z.  . (36) 

t t-1  t-z  t t-i 


of  Example  1. 
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Table  XXXI  shows  the  S and  R-arrays  using  f = (-l)'''p  for  a 

m m 

sample  of  size  25  generated  from  (36)  with  i.i.d.  N(0,1)  while 
Table  XXXII  show  the  calculated  values  of  D(n,m).  Clearly  from 
Table  XXXI  the  process  is  ARMA  (2,1).  Moreover  inspection  of 
D(n,m)  in  Table  XXXII  also  identifies  the  process  as  ARMA  (2,1), 
As  was  mentioned,  3 realizations  were  generated  at  each  of  the 
previously  listed  sample  sizes.  For  samples  of  size  25  in  this 
example,  one  realization  yielded  a much  more  distinct  ARMA  (2,1) 
pattern  than  those  of  Table  XXXI  and  one  yielded  a much  less 
distinct  pattern.  When  the  sample  size  was  as  large 

as  100  all  realizations  were  easily  identifiable. 

TABLE  XXXI 

rn  A 

S and  R Array  for  Example  11  with  f = (-1)  p 

m m 

for  a Sample  of  Size  25 


Sj^(m) 

S2(m) 

S^{m) 

m 

Rj^  (m) 

R2  (t") 

R3(m) 

-4.049 

4.667 

2.151 

-6 

-.082 

-.103 

-.032 

-2.383 

3.168 

-7.168 

-5 

.251 

.099 

-.053 

-1,704 

5.158 

84.758 

-4 

-.347 

-.085 

.021 

-1.112 

5.505 

5.958* 

-3 

.245 

.309 

.341 

11.343 

9.595* 

-1.235* 

-2 

-.028 

.459 

. 552* 

-3.936* 

1.558* 

2.896 

-1 

-.341 

-.659* 

.115 

-1.340* 

3.224 

-13.301 

0 

1.000* 

.107 

.012 

-.918 

3.084 

-1.097 

1 

-.341 

-.269 

-.039 

-9,867 

3.728 

4.  532 

2 

-.028 

-.185 

-.027 

-2,419 

2.714 

4.384 

3 

TABLE 

.245 

XXXII 

.100 

-.073 

D (n,m) 

From  a 

Sample  of 

Size 

25  in  Example 

11  at  f = 
m 

(-1)  Pj 

n 

0 

1 

2 

3 

4 

nT^ 

0 

8.391* 

.011 

.143 

.006 

.158 

1 

30.750* 

.000 

23.970 

.001 

1.852 

2 

.014* 

2.791 

.292 

,042 

.000 

3 

.712* 

.056 

.103 

.069 

.001 

The  symbol  * in  D(n,m)  tables  is  to  remind  the  reader  that 


n=0  is  not  included  in  the  D-Statistic  since  these  values  are  not 
scaled  the  same  as  other  values  of  n. 


•o  > 
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It  should  be  noted  that  the  pattern  in  the  S-array  is  much 
more  clear  if  one  focuses  on  the  bottom  part  of  S^Cm).  From  Monte 
Carlo  studies  made  by  the  authors,  this  is  to  be  expected.  That 
is,  in  some  limited  investigations  the  authors  found  S (-p+1) , 

A A P 

S (-p+2)  ...,  S (-p+n)  to  have  lower  mean  square  errors  than  their 
P A P * ^ 

counterparts  S (-p) , S (-p-1)  . . . , S (-p+n-1) , and  hence  for  very 
P P P 

small  samples  it  is  probably  wise  to  concentrate  on  the  region 
very  close  and  below  the  starred  values.  In  this  example  the  D- 
Statistic  clearly  gives  the  proper  estimate  for  p and  q.  Again, 
however,  for  samples  this  small  it  is  probably  better  to  use  a 
table  such  as  Table  XXXII  as  simply  a guide  for  visual  inspection. 
Finally  we  should  mention  that  neither  this  example  nor  the  ones 
that  follow  are  meant  to  suggest  that  25  samples  are  usually 
enough  data  to  determine  p and  q in  an  ARMA  (p,q)  process.  They 
do,  however,  suggest  that  in  many  cases  p and  q can  be  accurately 
determined  even  for  complex  processes  with  a surprisingly  small 
amount  of  data. 

Example  12. 

■ Consider  now  the  process  of  Example  6,  i.e. 

I (1-1.2B+-45B  ) (1-1,7B+.99B  )X  = (l-.SB)Z  (37) 

[ The  S and  R arrays  are  given  in  Table  XXXIII  for  a sample  of 

j size  25.  The  process  is  indicated  there  as  an  ARMA  (2,1).  More- 

over  it  is  clearly  nearly  nonstationary.  Proceeding  as  before 
I one  obtains  the  operator  1-1.69B+B^  as  the  proper  transformation 

to  stationarity . Upon  applying  this  operator  to  the  data  to  obtain 
= X^  - 1.69X^_j^  + only  the  moving  average  term  remained 

identifiable,  i.e.  the  process  appeared  as  an  MA  (1).  This  is  not 
surprising  since,  as  a review  of  Example  6 shows,  the  nonstation- 
ary factor  is  by  far  the  more  dominant  feature  of  this  model.  Al- 
though one  of  the  three  realizations  of  length  50  clearly  yielded 
the  proper  nonstationary  AP^  (4,1)  model,  a sample  of  size  200 
was  required  before  this  total  pattern,  i.e.  both  factors  in  (37), 
appeared  with  any  consistency. 


TABLE  XXXIII 


S and  R Arrays  for  a Sample  of  Size  25 


for  Example  12  at  f = (-D^'p 
m m 


Sj^(m) 

S2(m) 

S3(m) 

m 

Rj^(m) 

R2  (ni) 

R3  (m) 

-1.999 

3.845 

-13.475 

-6 

-.772 

-.125 

.000 

-1.675 

3.841 

5.553 

-5 

.771 

.162 

.000 

-1.154 

3.838 

39.622 

-4 

-.520 

-.377 

.001 

4.294 

3.833 

7.697* 

-3 

.080 

-.714 

.009 

-2.960 

4.004* 

-2.341* 

-2 

.423 

.211 

.027* 

-2.207* 

3.365* 

-3.233 

-1 

-.828 

-.171* 

.008 

-1.828* 

3.559 

-2.097 

0 

1.000* 

.144 

.001 

-1.510 

3.546 

-4.960 

1 

-.828 

-.196 

.000 

-.811 

3.547 

-35.552 

2 

.423 

.660 

.000 

-7.513 

3.546 

-7.302 

3 

.080 

.348 

.001 

Example 

13. 

In 

this  example 

we  consider 

a sample 

of  size  25 

from  the 

process 

of  Example 

9,  i.e. 

(l-ip 

^B)  (l-i|^2®> 

B)  (1- 

>1-48)  = Z^, 

(38) 

where 


= -.99,  rp2  ~ ’('3  ~ .5+41,  41^  = .5-.4i. 

The  S-array  at  f = (-l)"'p  is  shown  in  Table  XXXIV.  The  small 
mm  2 

values  of  S (m)  are  characteristic  of  the  factor  1-B  . Moreover 
in'* 

when  f = (-1)  p the  small  tabular  values  in  Column  3 strongly 
m m 

suggests  -1  is  a root  of  the  characteristic  equation.  Performing 
2 

the  operator  1-B  in  the  sequence  (1+B) (1-B)  yields  the  S-array 

at  f = p in  Table  X;CXV  for  AX^  = “ X,  , and  the  S-array  for 

m m ^ t t t-1 

= X..  - X..  - at  f = (-D^’p  in  Table  XXXVI.  Note  that  Table 
t c ni  n\  2 

XXXV  gives  strong  confirmation  of  the  diagnosis  that  1-B  is  a 

'V 

factor  of  (B) , i.e.  for  AX^,  s,  (m)  -2  at  f = p . But  Table 

t 1 mm 

XXXVI  then  identifies  the  resulting  data  as  AR  (2,0).  This  can 
also  be  seen  from  the  table  of  values  of  D(n,m)  in  Table  XXXVII. 

A plot  of  the  true  autocorrelation  and  the  estimated  autocorrela- 
tion is  shown  in  Figure  12  out  to  seven  lags  which  is  the  maximum 
number  required  to  calculate  the  S2 (m)  column  of  Table  XXXIV.  It 
is  interesting  to  note  that  although  the  estimates  are  poor  their 
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general  pattern  remains  that  of  the  true  autocorrelation  and 
hence  explains  why  p and  q are  clearly  distinguishable  in  this 
example. 


TABLE  XXXIV 

S-Array  at  f = p for  a Sample  of  Size  25 
m m 


for  Example  13 


m 

Sj^(m) 

$2  (m) 

S^Cm) 

-6 

-.393 

- 

.150 

.639 

-5 

1.207 

- 

.115 

.448 

-4 

-.742 

- 

.222 

.279 

-3 

4.628 

- 

.001 

.279* 

-2 

-1.036 

- 

.278* 

- 

.122* 

-1 

-36.451* 

.219* 

- 

.146 

0 

1.028* 

.001 

- 

.145 

1 

-28.919 

.145 

- 

.112 

2 

-.822 

.073 

- 

.210 

3 

2.880 

.092 

- 

.308 

4 

-.547 

.055 

- 

.262 

TABLE  XXXV 

S-Array  for  a 

Sample  of 

Size 

25  for  AX 

t 

in  Example  13 

at  f 
m 

= Pm 

m 

S^(m) 

S2(m) 

-6 

-2.093 

-.490 

-5 

-2.023 

.029 

-4 

-2.105 

2.226 

-3 

-2.029 

-33.049 

-2 

-2.022 

-3.001* 

-1 

-2.139* 

1.156* 

0 

-1.878* 

1.858 

1 

-1.978 

21.525 

2 

-1.972 

-.027 

3 

-1.905 

.373 

64 


TABLE  XXXVl 

S-Array  for  a Sample  of  Size  25  For 

W = X -X  _ in  Example  13  at  f = (-l)"'p 
t c m r 


m 

Sj^(m) 

S2(m) 

SjCm) 

(m) 

-6 

-2.309 

-4.028 

-9.327 

13.404 

-5 

-2.799 

6.888 

1.022 

13.036 

-4 

-1.993 

4.197 

.386 

13.260* 

-3 

-.572 

4.519 

9.343* 

4,267* 

-2 

-6.613 

5.707* 

-1,792* 

4.968 

-1 

-2.675* 

2.218* 

-.183 

4.979 

0 

-1.597* 

2.887 

-.429 

.469 

1 

-1.178 

2.324 

-3,981 

.628 

2 

1.338 

2.189 

-.203 

10.378 

3 

-2.007 

1.041 

-2.462 

-10.649 

4 

-1.556 

-.287 

-3.029 

-8.828 

TABLE  XXXVII 

Values  of 

D(n,m)  in  Example 

13  For 

”t  = \ 

-X  _ at  f = 

(-1) 

m* 

% 

\n 

0 1 

2 

3 

4 

0 

2.201*  .017 

1.970 

.002 

.001 

1 

2.180*  .003 

.008 

.356 

.000 

2 

.099*  .779 

.002 

.001 

2.027 

As  is  clear  from  the  above 

tables  an  ARI-JA  (4, 

2)  is  also 

a possible 

choice. 

However  the  ARMA  (2,0)  pattern 

does 

loo)c  better 

in  Table 

XXXVI.  Moreover  with  a sample  this  small  one  would  undoubtly 
choose  the  ARMA  (2,0). 
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VIII.  SEASONAL  DATA 

Thus  far  we  have  not  considered  any  examples  where  there  was 
a strong  seasonal  component  of  rather  long  period  such  as  in  a 
model  of  the  type 

(l-B^^)X^  = (1-0B^^)2^  (38) 

At  first  glance  one  might  not  expect  the  methods  of  this 
paper  to  be  of  particular  use  in  this  problem,  due  to  the  large 
values  of  p and  q,  unless  a large  data  set  is  available.  The 
extent  to  which  this  is  true  is  not  presently  known  to  the  authors. 
However,  a number  of  modifications  can  be  made  to  the  methods  of 
this  paper  which  might  allow  a general  approach  to  the  seasonal 
problem.  In  case  of  (38)  such  modifications  are  easily  found  so 
that  handling  that  particular  case  is  of  no  difficulty.  However 
the  authors  have  yet  to  establish  a general  approach  to  the  prob- 
lem at  this  time.  Some  preliminary  results  have  been  established 
and  can  be  found  in  Kelley  (1977)  where  an  interesting  example  of 
the  so  called  "Colored  Fox  Data"  (see  Anderson,  1977)  is  given. 
Hopefully  a comprehensive  approach  to  the  problem  will  appear  in 
a later  paper. 

IX.  THE  MOVING  AVERAGE  PROCESS 

Thus  far  only  cases  where  p > 0 have  been  considered.  This 
is  not  to  suggest  that  the  method  is  not  useful  if  the  data 
happens  to  come  from  an  MA  (q)  process.  In  fact  these  authors 
have  found  that  to  be  one  of  simpler  cases  to  distinguish.  How- 
ever, one  should  be  reminded  that  in  the  case  of  a MA  (q)  the 
S-array  will  have  no  fixed  pattern  and  hence  the  method  will 
direct  one  back  to  the  R-array  which  in  turn  leads  one  to  estimate 
q from  R^(m),  i.e.  from  the  values  of  the  autocorrelation.  Thus 
at  this  point  the  method  coincides  with  that  of  Box  and  Jenkins. 
The  D(n,m)  values  at  n = 0 are  useful  in  aidina  one  in  not  over- 
looking the  MA  (q)  possibility  when  inspectina  the  R and  S arrays- 

However,  only  comparing  them  to  other  D(n,ml  values  would  lead  to 
the  choice  of  a moving  average  too  often. 
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X.  NON-NORMAL  NOISE 

In  all  previous  examples  was  i.i.d.  N(0,1).  Of  course 
since  the  statistic  of  interest  was  p there  was  no  loss  in  gen- 
erality  in  using  only  unit  variance.  However,  the  question  of 
how  much  the  method  would  be  effected  if  the  noise  were  not 
normal  is  of  interest.  Thus  the  authors  considered  all  of  the 
previous  examples  again  for  samples  of  the  same  size  as  initially 
considered  but  with  first  exponential  noise  and  then  Cauchy  noise. 
For  these  examples  and  saimple  sizes  no  detectadjle  difference  was 
found  in  the  effectiveness  of  the  method.  For  small  samples  this 
may  not  be  the  case  but  that  was  not  investigated. 

XI.  CONCLUDING  REMARKS 

In  this  paper  the  authors  have  presented  what  they  feel  is  ^ 
satisfactory  approach  to  the  problem  of  estimating  p and  q in 
ARMA  (p,  q)  modeling  problems.  Several  additional  topics  con- 
cerning this  problem  have  been  investigated  but  not  considered 
here,  the  most  promising  being  seasonal  models  and  the  transfer 
function  model.  Work  in  these  areas  will  appear  in  later  papers. 

Some  closing  remarks  should  be  made  concerning  the  D- 
Statistic.  The  current  definition  may  be  a bit  ad  hoc  and  some 
modifications  of  the  statistic  may  be  necessary.  At  the  present 
this  was  as  good  as  the  authors  could  do  in  the  way  of  finding 
a single  measure  of  the  presence  of  the  desired  pattern.  Un- 
doubtly  some  improvement  can  be  made  there.  In  addition  it 
seems  that  something  similar  to  the  D-Statistic  could  be  done 
in  the  nonstationary  case  so  that  the  data  could  be  transformed 
to  stationarity  "rather  automatically"  if  necessary. 

Finally  it  should  be  mentioned  that  programs  for  calculating 
R,  S and  D(n,m)  arrays  are  included  in  the  Appendix  so  that 
application  of  the  methods  of  this  paper  should  not  be  difficult. 
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APPENDIX  1 

Proof  of  Theorem  3:  Since  is  stationary  ARMA  (p,  q) 

e L{p,A)  for  p > q + 1.  Then  there  exist  constants 

Ip. , , p such  that 
Ip  p 


p = E p.p.  . for  k > q + 1. 
^ i=l  ^ 


Moreover 


1 

1 

1 

p-k 

P-k+1  • • • 

P-k+p 

•^-k+l 

P-k+2  ' ' • 

P-k+p+1 

^-k+p-1 

P-k+p  • • ' 

P-k+2p-l 

P-k 

P-k+1  ■ ■ ‘ 

P-k+p-1 

P-k+1 

P-k+2  • • • 

P-k+p 

P-k+p-1 

0 , , • • • 

-k+p 

P-k+2p-2 

is  not 

zero. 

In  light  of  (Al)  and  the  fact  that  the  autocorrelation  function 
is  an  even  function  of  the  lag,  the  numerator  of  (A2)  may  be  written 

1 P P 

- ^ (-1)^(1-  .1.  P.)  H [p  .]  . 

<tip  1=1  1 p -k 

for  -k  < -p  -q,  by  applying  appropriate  row  and  column  operations. 
Then  o , v ^ ^ 

S (P.v)  = - r~  ^ 

p -k  p^ 

and  C 5^  0 since  the  process  is  stationary. 

Suppose  n is  the  smallest  positive  integer  such  that 

^ “ 


for  -k  < M.  Then  from  (a2) 


1 1 ...  1 

0 0 ...  0 1 

P-k  ^-k+1  . . . P-k+n 

• • 

• • 

• • 

= (-d'^C 

P-k  P-k+r  " P-k+n-1  P-k+n 

« • • • 

• • • • 

. • • • 

P-k+n-1  P-k+n  ‘ ' ' P-k+2n-l 

P-k+n-1  P-k+n"  •P-k+2n-2  P-k+2n-l 
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so  that 


1 

1 

1 

i-(-i)''c 

p 

P 

• • • 

P 

P 

-)c 

-k+1 

-k+n-1 

-k+n 

• 

« 

• 

’ 

• 

'^-k+n-l 

^-k+1 

• • • 

’^-k+2n-2 

'’-k+2n-l 

= 0. (A4) 


That  is  a determinant  such  as  (A4)  may  be  exhibited  for  each  k 

such  that  -k  < M.  Then  using  the  symmetry  of  p , for  each  k 

constants  b, (k) , b- (k) , . . . , b^ (k)  exist  such  that 
1 z p 


■ 1 

1 

■i-(-i)"c* 

tr 

^k 

+ ..  . + b (k) 
n 

‘’k-n+l 

+ 

^k-n 

-^k-n+1  • 

. ^k-2n+2. 

.‘^k-2n+l  . 

0. 


Consequently 


■ 1 1 ...  1 ■ 

■b^(k)-bj^(k+l)' 

^k  ‘’k-l  •••  ‘^k-n+l 

• • • 

• 

• • • 

• • • 

0*  0*  4 ••• 

b (k)-b  (k+1) 

‘^k-n  ^k-n-1  k-2n+l 

n n 

^ • 

. 

= 0.  (A5) 


Then  either  the  determinant  of  the  nxn  matrix  of  (A5)  is  zero  or 

the  vector  is  null.  However  the  determinant  of  the  nxn  matrix  is 

the  numerator  of  S , (p  ) . That  this  quantity  would  be  zero  is  a 
n-1  m 

contradiction  to  the  assumption  that  n is  the  smallest  value  such 
that  (A3)  holds. 

Hence 

b.  (k)  = b.  (k+1)  for  i = 1,  2,  ...,  n 
1 1 

amd  k > -M  . 

Thus  the  constants  b^  ( k ) do  not  depend  on  k and  we  have 

n 

^k-n+1  ’^i*^k-n+l-i 

1=1 

for  k-n+1  > -M-n+1.  It  follows  therefore  that  n = p and  M = -q-p, 


completing  the  proof. 
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APPENDIX  2 


(x,y) : X = minimum  lag  of  acf,  y = maximum  lag  of  acf 
required  to  compute  n=l,2,.,.,  lc=. . . ,-l,0,l,2,  ... 


k " 

1 

2 

3 

4 

5 

-12 

(11,12) 

(9,12) 

(7,12) 

(5,12) 

(3,12) 

-11 

(10,11) 

(8,11) 

(6,11) 

(4,11) 

(2,11) 

-10 

(9,10) 

(7,10) 

(5,10) 

(3,10) 

(1,10) 

-9 

(8,9) 

(6,9) 

(4,9) 

(2,9) 

(0,9) 

-8 

(7,8) 

(5,8) 

(3,8) 

(1,8) 

(0,8) 

-7 

(6,7) 

(4,7) 

(2,7) 

(0,7) 

(0,7) 

-6 

(5,6) 

(3,6) 

(1,6) 

(0,6) 

(0,6) 

-5 

(4,5) 

(2,5) 

(0,5) 

(0,5) 

(0,5)* 

-4 

f3,4) 

(1,4) 

(0,4) 

(0,4)* 

(0,5)* 

-■) 

(2,3) 

(0,3) 

(0,3)* 

(0,4)* 

(0,6) 

-2 

(1,2) 

(0,2)* 

(0,3)* 

(0,5) 

(0,7) 

-1 

(0,1)* 

(0,2)* 

(0,4) 

(0,6) 

(0,8) 

0 

(0,1)* 

(0,3) 

(0,5) 

(0,7) 

(0,9) 

1 

(1,2) 

(1,4) 

(1,6) 

(1,8) 

(1,10) 

2 

(2,3) 

(2,5) 

(2,7) 

(2,9) 

(2,11) 

3 

(3,4) 

(3,6) 

(3,8) 

(3,10) 

(3,12) 

4 

(4,5) 

(4,7) 

(4,9) 

(4,11) 

(4,13) 

5 

(5,6) 

(5,8) 

(5,10) 

(5,12) 

(5,14) 

i 
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(x,y) : X * minimum  lag  of  acf , y = maximum  lag  of  acf 
required  to  compute  r^(k),  n = 1,2,...,  k = ...,-1,0,1,... 


k" 

1 

2 

3 

4 

5 

-12 

(12,12) 

(10,12) 

(8,12) 

(6,12) 

(4,12) 

-11 

(11,11) 

(9,11) 

(7,11) 

(5,11) 

(3,11) 

-10 

(10,10) 

(8,10) 

(6,10) 

(4,10) 

(2,10) 

-9 

(9,9) 

(7,9) 

(5,9) 

(3,9) 

(1,9) 

-8 

(8,8) 

(6,8) 

(4,8) 

(2,8) 

(0,8) 

-7 

(7,7) 

(5,7) 

(3,7) 

(1,7) 

(0,7) 

-6 

(6,6) 

(4,6) 

(2,6) 

(0,6) 

(0,6) 

-5 

(5,5) 

(3,5) 

(1,5) 

(0,5) 

(0,5) 

-4 

(4,4) 

(2,4) 

(0,4) 

(0,4) 

(0,4)* 

-3 

(3,3) 

(1,3) 

(0,3) 

(0,3)* 

(0,5) 

-2 

(2,2) 

(0,2) 

(0,2)* 

' (0,4) 

(0,6) 

-1 

(1,1) 

(0,1)* 

(0,3) 

(0,5) 

(0,7) 

0 

(0,0)* 

(0,2) 

(0,4) 

(0,6) 

(0,8) 

1 

(1,1) 

(1,3) 

(1,5) 

(1,7) 

(1,9) 

2 

(2,2) 

(2,4) 

(2,6) 

(2,8) 

(2,10) 

3 

(3,3) 

(3,5) 

(3,7) 

(3,9) 

(3,11) 

4 

(4,4) 

(4,6) 

(4,8) 

(4,10) 

(4,12) 

5 

(5,5) 

(5,7) 

(5,9) 

(5,11) 

(5,13) 
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APPENDIX  3 


Theorem  2A;  Let  and  be  the  roots  of  the  characteristic 

equation  1 - (j>^x  - <^2^^  ~ = 0 of  an  ARMA  (3,1)  process  in 

which  M,  1 < M < 3,  of  the  roots  r^,  r2  and  r^  lie  on  the  unit 

circle.  Let  the  integer  k,  1 < k < M be  the  number  of  distinct 

roots  having  modulus  1.  Then  p*  as  defined  in  Definition  6 satis- 

m 

fies  a linear  homogeneous  difference  equation  with  constant  co- 
efficients whose  order  is  less  than  or  equal  to  k for  m = 0, 

±1,  and  whose  characteristic  equation  has  all  of  its  roots 

distinct  cuid  on  the  unit  circle. 

Preliminaries  to  Proof;  For  a stationary  ARMA  (3,1)  process, 


Now 

Pm  = '^’iPm  1 + ? + 

m 1 m— X X m—x 

the  general  solution  to  (1)  is 

4),P  , for  m > 2. 

j m- j - 

(1) 

p = A + A r~l"*l  + A r~ 

m 1 1 2 2 3 3 

for  3 distinct  roots. 

(2) 

and 

"m  = (^1  ^ ^ 

^ for  r^^  repeated  once 

(3) 

•^m  = '‘=1  ^ S'"  ^ 

for  r^^  repeated  twice. 

(4) 

In  these  equations  A^,  A^,  A^,  B^,  B^,  B^,  C^,  C2 , and  are 
all  constants  determined  by 


where  letting  A = 
B = 


I and 


1 

l-<ji2~'i'3  (<(>]_+t3)  1*^1  ' 

<(>1  + + (<f'2+'l'3ii>3)  (<(>3+<t>3)  and 

1 ~ 4^2  “ 'i* 3 ^ 3)  ' 

(<i>^+((>243)A  - [l-<!i3-42(4>2+<(>3_<(!3)  IB 
B{B[l-6^(<).^-0^)]-9j^A}  ' 


I 


where 


and 
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^ = (G2-G3)(p^(G2+G3)  - P2  - PoG2G3)/D 


(G3-G^) (P3(G3+G^) 


^2  " 


B, 

X 


(G^-G^Xp^CG^+G^)  - P2  - 

(Gi-G^) (G^-G3) (G3-G3)  ; 

G3  ^ - P2 

(Gi-G3)^ 


P2  - P3^(G^+G3)  + Gj^G3 


PoV2’/° 


Gi(G^-G3) 


P2  - 2G3P3_  + Gj_ 


(G^-G3)' 


<=1=  ^ 


C,  = 


4P1G3 


P2  - 3G^ 


2 

- 2Pj^G3  + G3 


Finally  recall 

^1  = Gi  + G2  + G3 

4-2  = - (G^G2  + G^G3  + G2G3 

(t>,  = G G G 
^3  12  3 

Essentially,  the  proof  involves  expressing  the  appropriate 
general  solution  (2) , (3) , or  (4)  of  the  difference  equation  (1) 
in  terms  of  G^^,  G2  and  G3  and  taking  the  limit  of  the  general 
solution  as  the  specified  roots  approach  the  unit  circle.  In  each 
case,  the  roots  which  approach  the  unit  circle  are  arbitrarily 
chosen.  We  have  the  following  results: 
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Proof  of  Theorem  2A: 


Case  1;  M = 1,  k = 1,  all  roots  real.  Let  -►  1 while 

and  G remain  fixed  with  |g,|  <1  and  |g,|  < 1.  Then  lim  A = ] 

^ ^ G^-1  ^ 


lim  A 


2=0  and  lim  A - 0.  Also,  lim  p = 1 and  lim  p = 1. 


Thus,  the  acf  of  the  process  is  given  by  p*  - 1 for  m ■ 0,  ±1, 

m 

. . . ; i.e.,  p*  = p*  , and  we  see  that  p*  is  a solution  to  a first 
m m— 1 m 

order  difference  equation. 

Case  2:  M = 2,  k = 1,  all  roots  real.  Assume  G^^  has  multi- 


plicity 2 and  take  the  limits  as  G, 


1 where  iGjI  < 1 with  G^ 


fixed.  The  results  are  identical  to  the  results  given  in  Case  1; 


i.e.,  p* 
m 


p*  , for  in  = 0,  ±1, 
in— ± 


Case  3:  M = 3,  k = 1,  all  roots  real.  Assuming  G^^  has 
multiplicity  3,  these  results  are  also  identical  to  those  in  Case 


Case  4:  M = 2,  k = 2,  2 complex  roots,  1 real  root.  Let 

i 0 —10 

= X + jy  = de"^  and  G2  = G^  = x - jy  = de  ^ where  d = 

/x  + y , j = and  tan  9 = x/y.  Assume  G^  is  fixed  with 

|g,  I < 1.  The  limits  are  taken  as  d 1 . Then,  lim  A.  = 

d-*-r  ^ 

lim  A = 1/2,  lim  A = 0,  lim  p^  = cos  9,  and  lim  P,  = 

d-*-l“  d-*-l“  d-*-l“  d-<-l“ 


cos  29.  Thus,  the  general  solution  for  the  difference  equation 

, .1  j9|m|  , 1 -j9|m|  q „ m ^ 

becomes  p*  = -e  ' - e ' ' = cos  9m  for  in  = 0,  ±1,  ...  , 

m2  2 

which  is  the  solution  to  the  second  order  difference  equation 

whose  characteristic  equation  has  roots  e""  and  e •'  . 

Case  5;  M = 2,  k = 2,  all  roots  real.  Here,  we  have  Gj^  -*•  1, 

G.  -*■  -1  with  G-  fixed  and  |g,1  < 1.  The  limits  are  taken  with 
2 3 'I  3'  ^ ^ 

the  restriction  G2  = -G^.  Then  lim  A^^  = -(p?+l),  lim  A2  = r(l-pj), 

G^-*-!  Gj^-*-l 

lim  A = 0,  where 

Gi^l 
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I 


p*  = lim  p = = 

l+G^ 


2 2 

20^(1-03)^ 


(1+G3)  t (i-e^  ^ ~S®i  ^ 


Also,  lim  p.  = 1. 

Gi-1 

equation  becomes 


Thus,  the  general  solution  to  the  difference 


P*  = 5t(P*+l)  + (-1)  (1-P*)] 

m 1 1 


m odd 


, m even 


i.e.,  p*  = p*  for  M = 0,  ±1, 
rti  nv— 4 


, which  implies  that  p* 

m 


satisfies  a second  order  difference  equation  with  roots  ±1. 

Case  6;  M = 3,  m = 2,  all  roots  real.  Let  G^^  have  multi 

The  limits  are  taken  with  the 


,plicity  2,  G^  1,  and  G^  -<■  -1. 


restriction  G^  = -G^^. 


Then 


lim  B. 


lim  B.  = 1,  lim  B. 

G^->1  Gj,-*-l  ^ Gj^-^-l 


0, 


lim  p, 

Gi^l  ^ 


lim  p = 1.  These  results  are  seen  to  be  identical  to 

Gi-1 

those  in  Case  1;  i.e.,  p*  = p*  , for  m = 0,  ±1,  ...  . (This  case 

m m“i 

is  distinct  in  that  the  order  of  the  resulting  difference  equation 
is  less  than  m,  the  number  of  distinct  roots  on  the  circle,  while 
in  the  preceding  cases,  the  order  was  equal  to  m. ) 

Case  7:  M=3,  m=3,  2 complex  roots,  1 real  root.  Let 

= de  ' 


G^  = de 


j9 


and  take  the  limits  as  G. 


1, 


and  G^  -»■  1. 


Under  the  restriction 
d = G,  and  take  the  limits  as  G, 


IGJ  = iG^i 

1 . Then : 


= G. 


1-2' 
we  let 


lim  A,  = lim  A,  ^ (rr-^ — S 
G3-I  ^ G3-I  ^ 2 2+cos  9 


, , • , 1+cos  9 

' 3 ~ 2+cos  9 

G3-I 


lira  P 

G3-I 


2 cos  9-H 
cos  9+2 


lira  p * 
G3.1 


cos  9 (2  cos  9+1) 
2+  cos  9 


Thus,  p*  is  a 
in 


solution  to  a third  order  difference  equation  for  ra  = 0+1,  +2  ...  . 

Other  cases  are  possible  but  have  been  omitted  here  because 
of  their  near  equivalence  to  one  of  the  previous  cases.  For  exam- 
ple if  M = 1,  k = 1 and  G^^  -*•  -1  the  result  is  the  same  as  Case  1 

except  p*  = -p* 

m m-1 
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SUflROUTIN'E  RANDS  (MXHOWfMROW.N'  G * KOL  * I PHO . G * R . S t I WK  . I R S ) 


DESCRIPTION  OF  PARAMETERS 

MXROW  - ROW  DIMENSION  OF  R AND  S ARRAYS  IN  CALLING  PROGRAM. 

(INPUT) 

MROW  - number  of  rows  OF  R AND  S ARRAYS  TO  BE  CALCULATED  (INPUT) 
NEG  - NUMBER  OF  ROWS  OF  NEGATIVE  LAG  IN  R AND  S ARRAYS 
NEG  < MROW  - 1 
( INPUT) 

KOL  - NUMBER  OF  COLUMNS  OF  P AND  S ARRAYS  TO  BE  CALCULATED 

R AND  S ARRAYS  MUST  BE  DIMENSIONED  TO  AT  LEAST  KOL  COLUMNS 
IN  CALLING  PROGRAM.  (INPUT) 

KOL  < (MROW  ♦ 2)/?  AMO  0 < KOL  < 12 

THE  RFSTRICTION  THAT  KOL  IS  LESS  THAN  TWELVE  IS  IMPOSED  BY 
THE  WRITE  STATEMFrjTS.  WITH  IRS  .NE.  0.  KOL  MAY  BE  AS 
LARGE  AS  PERMITTED  BY  MROW. 

IPHO  - IRHO  .EO.  0.  R AND  S CALCULATED  FOR  AUTOCORRELATIONS  IN 
COLUMN  1 OF  THE  R ARRAY, 

IRHO  .NE.  0.  R AND  S calculated  For  (-1)<k>LAG  TIMES  THE 
AUTOCORRELATION  FUNCTION  IN  COLUMN  1 OF  THE  R ARRAY. 

( INPUT) 

G - ARRAY  OF  AUTOCORRELATION  ESTIMATES  OF  LAG  ZERO  TO  N FOR 
THE  FIRST  COLUMN  OF  THE  R ARRAY.  WHERE  N IS  AT  LEAST  THE 
MAXIMUM  OF  MROW-NEG-I  AND  NEG  . (INPUT) 

R - R ARRAY.  DIMENSIONED  MXROW  BY  K IN  CALLING  PPOGRAMt 
where  K is  greater  than  or  equal  to  KOL.  (OUTPUT) 

S - S ARRAY.  dimensioned  AS  R ARRAY  IN  CALLING  PROGRAM. 
(OUTPUT) 

IWK  - WORK  VECTOR  DIMENSIONED  AT  LEAST  KOL  IN  CALLING  PROGRAM. 

IRS  - PRINT  CONTROL  PARAMETER.  (INPUT) 

IRS  .EO.  0i  R AMO  S ARRAYS  ARE  PRINTED  ON  UNIT  NUMBER  IPR» 
WHICH  IS  SET  IN  the  data  STATEMENT. 

IRS  .NE.  0,  CAUSES  PRINTING  TO  BE  SUPPRESSED. 


DIMENSION  R (MXpnW. 1 ) .S  (MXRDW, 1 ) t G ( 1 ) . I WK  ( 1 ) 
data  IBLK/IH  /, ISTR/IH*/. IPR/6/, IZRO/0/ 
KOLMl=KOL-' 

DO  5 IC=l.KOL 
no  5 IR=l.MROW 
R(IR.IC)=0.0 
5 S(IR. 10=0.0 
MN=MP0W-NEG 
DO  10  I=1.MN 
NEGI=NEG*I 
S (NEGI . 1 ) =1 .0 
10  R (MEGI » 1 ) =G ( I ) 

DO  20  I=1.NFG 
NEGI=NEG*2-I 
S ( 1 . 1 ) =1  . 0 
P(I»1)=G(NERI) 
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IFdPHO  .EO.  0)  GO  TO  40  I 

DO  30  I=ltMROW  I 

NEGI=r-NEG-l  j 

30  R(I»l)=R(Itl)*(-U**(N'EGI)  : 

40  Jl=MROW 

DO  70  IC=2.K0L 

ICM1=IC-1 

JlsJl-1 

DO  60  IR=ltJl 

IRP1=IR*1 

IF ( ABS (R ( IR» ICMl ) ) .EQ.  0,0)  GO  TO  50 
S{IR,IC)=S(IPP1»ICM1)«(R(IRP1,ICM1)/R(IR*ICm1)-1,0) 

GO  TO  60 

50  S(IR.IC)=-S(IRPltICMl) 

60  CONTINUE 
J1=J1-1 
DO  70  IR=1,J1 
IRP1=IR*1 

70  R(IR.IC)=R(IRPldCMl)  »(S(IRP1»IC)/S(IR*IC)-1.0) 

IF (IRS  .NE.  0)  RETURN 
IFdRHO  .EO.  0)  GO  TO  80 
WRITE  dPR»200) 

GO  TO  90 

80  WRITE  dPR. 210) 

90  write  dPP»220)  dCdC  = l«KOL) 

DO  110  IR=l,PROW 
IPN=IR-NEG-1 
DO  100  IC=l.KOL 
IWK ( IC) =I8LK 

100  IFdR  ,E0.  NEG*2-IC)  IWKdC)=ISTR 

no  WRITE  dPR^230)  IRNi  ( R (I  R , K ) d WK  ( K ) » K = I » KOL ) 

WRITE  dPRi 24  0 ) IZROi  (ICi IC=1 .KOLMl ) 

DO  130  IP=1,nROW 
IRM=1P-NEG«1 
DO  120  IC=nKOL 
IWK ( IC) =I8LK 

IFdR  .EO.  NF6+2-IC)  rWKdC)=ISTR 
120  IFdR  .EO.  NEG  + 3-IC)  lWKdC)=ISTR 
130  WRITE  dPR. 230)  IRN,  ( S (I  R » K ) * I WK  ( K ) » K=  1 1 KOL) 

RETURN 

200  FOPMAT(//20X.27H(  F(M)  = ( ( - 1 ) ) « ACF  ( M ) )) 

210  F0RNAT(//25X,17H(  F(M)  = ACF(M)  )) 

220  FORMAT (//30X ,7hr-ARRAY//3H  K.11(9X.I2)) 

230  F0RMAT(2H  (,I3.1H)  .11  dX,F9.4,Al)  ) 


240  FOP^^AT  (//30x,7HS-ARRAY//3H  K,11(9X.I2)) 

END 
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SUBROUTINE  DSTAT ( MXHOW , MROW .MEG i R . S * IP. IQt IPNT) 


DESCRIPTION  OF  PARAMETERS 

MXROW  - ROW  DIMENSION  OF  R AND  S ARRAYS  IN  CALLING  PROGRAM, 

(INPUT) 

MROW  - NUMBER  OF  ROWS  TO  WHICH  R AND  S WERE  CALCULATED.  (INPUT) 
NEG  - NUMBER  OF  ROWS  OF  NEGATIVE  LAG  IN  R AND  S ARRAYS 
(INPUT) 

R - R ARRAY,  DIMENSIONED  MXROW  BY  K IN  CALLING  PROGRAM. 

WHERE  K IS  GREATER  THAN  OP  EQUAL  TO  KOL.  (INPUT) 

S - S ARRAY.  DIMENSIONED  AS  R ARRAY  IN  CALLING  PROGRAM, 

( INPUT) 

IP  - (INPUT)  MAXIMUM  ORDER  OF  AUTOREGRESSION  PERMITTED. 

(OUTPUT)  ORDER  OF  AUTOREGRESSION  SELECTED. 

IQ  - (INPUT)  MAXIMUM  ORDER  OF  MOVING  AVERAGE  PERMITTED. 

(OUTPUT)  ORDER  OF  MOVING  AVERAGE  SELECTED. 

IPNT  - PRINT  CONTROL  PAREMETER  (INPUT) 

IPNT  .EO.  0,  D statistic  VALUES  ARE  PRINTED, 

IPNT  .ME.  0.  PRINTING  IS  SUPPRESSED. 

WITH  IPNT  = 0 THE  WRITE  STATEMENTS  REQUIRE  THE  INPUT  VALUE 
OF  10  TO  BE  less  THAN  TEN. 

PRINTING  IS  ON  UNIT  IPR  WHICH  IS  SET  IN  THE  DATA  STATEMENT. 

NOTE 

NEG  MUST  BE  GREATER  THAN  OR  EQUAL  TO  IP*IQ*3.  FOR  IP  AND  IQ 
THEIR  INPUT  VALUE‘S.  FURTHERMORE  MROW  MUST  BE  GREATER  THAN  OR  EQUAL 
TO  NEG.IP*IQ*a,  for  ip  and  IO  EQUAL  TO  THEIR  INPUT  VALUES. 

IF  ( IP*  1 ) *»  ( IQ.  1 ) IS  GREATER  THAN  SO  (WHERE  IP  AND  IQ  ARE  AT 
THEIR  INPUT  VALUES)  THEN  THE  DIMENSION  OF  STAT  MUST  BE  CHANGED 
ACCORDINGLY. 


DIMENSION  R (MXROW,  1 ) ,S (MXROW,  1 ) ,STAT  (50 ) 
DATA  IPR/6/, I ZRO/0/ 

JPP=IP, 1 
JOQ=IQ.  1 
rOP=JPP*JOQ 
DO  AO  1=1, JPP 
JP=I-1 

DO  AO  J=l.JQO 
JO=J-l 

JK  = J(DO*  ( I-l  ) .J 

ISUB2=NE6-JP.JQ.2 

ISU81=ISUB2-1 

IPl=I.l 

ru  = i-i 

ROT= (S ( ISUB2. I ) *5 ( ISUBl . IPl ) ) o»2 
IF(I  .GT.  1)  GO  TO  20 
DO  10  JI=1,3 
NEGM=NEG-JP-JQ-J1 ♦ 1 
NEGP=NEG-JP. JQ. JI  ♦ 1 
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10  80T  = R0T*f?  (NEOM,  I ) «*2  ♦ R{NEGP»I)**2 
GO  TO  35 

?o  isun3=isue2*i 

nOT  = BOT/(S(ISUB3,IMl)  ♦ S ( I SU82 f I ) ) *«2 
DO  30  JI=1»3 
NEGM=NEG-JP-JO-JI ♦! 

MEGP=NEG-JP*JO*jr  + I 
NEGMl =NEGM* I 
MEGP1=NEGP*1 

30  POT=BOT* (R (NEGM, I ) /R ( NEGM 1 , I M 1 ) ) (R (NEGP» I ) /R (NEGPl . IMI ) ) **2 
35  NEG0=NEG-JP-JQ 
HEG1=NEG0*1 

ARG  = S <ME60»IPn  / (80T*S(NEG1 , 1 ) ) 

AO  STAT ( JK) =A8S (ARG) 
iyx  = i 

XMX  = STAT ( 1 ) 

00  50  1=1 . lOP 

IF(STAT(I)  .LT.  XMX)  GO  TO  50 
IMX=I 

X.mx  = STAT  ( I ) 

50  CONTINUE 

IP=(IMX-l)/JOO 
IQ=lMX-(rP»JQQ)-l 
IFdPNT  .NE.  0)  return 
WRITFdPR,  100) 

JQO=JOO-l 

WRITFdPRtllO)  IZRO*  (I  d = l ♦ JQQ) 

DO  GO  J=1»JPP 
JM1=J-1 

IST=JM1*(J0Q>1) ♦! 

IEND=IST*JQQ 

GO  V/PITE  dPPd20)  JMl  . (STAT  (I ) , I = IST.  lENO) 

WRITE (IPR* 130)  IP, IQ 
RETURN 

100  FOPMAT(///50X,11HO  STATISTIC//2qX,llH0R0ER  OF  MA) 

110  FOPMATdGH  ORDER  OF  AR  , 1 0 ( 4 X , T 3 , 4 X ) ) 

120  FORVAT  (/4X  , r2»6X,10dX,E10.4)) 

130  F0RMAT(/34H  order  OF  AUTOREGRESSION  SELECTED , 2X , I 4/ 

1 34H  ORDER  OF  MOVING  AVERAGE  SELECTED  , 2X , 1 4 ) 

END 
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